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ABSTRACT 


A  six  degree  of  freedom,  20-state  model  of  two  spacecraft  rendezvous  is 
developed,  one  of  which  was  controlled  and  the  other  considered  to  be  passively 
tumbling.  Solutions  that  minimize  a  series  of  performance  indices  are  obtained  for  the 
problem  of  close  approach,  up  to  the  point  of  contact,  using  a  direct  optimal  control 
method.  The  solution  is  then  verified  as  optimal  by  way  of  an  indirect  method  based  on 
the  Minimum  Principle.  Next,  a  trajectory  generation  method  for  spacecraft  reorientation 
is  developed,  based  on  a  quaternion  construction  of  Inverse  Dynamics  in  the  Virtual 
Domain.  This  new  construction  enables  development  of  an  Inverse  Dynamics  in  the 
Virtual  Domain  rapid-trajectory  generation  method  that  exploits  the  concept  of 
decoupling  space  and  time,  for  the  problem  of  a  spacecraft  perfonning  a  close  approach 
maneuver  to  a  tumbling  object.  Finally,  the  advantages  of  the  new  method  are 
demonstrated  through  simulated  scenarios  that  employ  two  distinct  concepts  of  closed- 
loop  feedback.  The  benefits  seen  by  Inverse  Dynamics  in  the  Virtual  Domain  methods 
include  the  rapid  computational  time  that  allows  a  feasible  solution  to  be  generated, 
potentially  onboard  a  spacecraft  and  in  closed-loop.  Although  the  Inverse  Dynamics  in 
the  Virtual  Domain  method  cannot  match  the  true  optimal  solution,  it  has  several 
advantages.  Rapid  computational  time  and  the  ability  to  reshape  itself  when  provided 
updated  state  variable  information  reinforce  the  overall  robustness  of  the  method  for  safe 
trajectory  planning.  The  novel  trajectory  generation  method  developed  is  tested  using 
Monte  Carlo  methods  to  demonstrate  its  ability  to  handle  realistic  situations  with  varying 
initial  conditions. 
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I.  INTRODUCTION 


The  rendezvous  problem  of  two  spacecraft  orbiting  the  earth  has  been  addressed 
in  numerous  publications,  research  thrusts  and  desired  mission  capabilities.  Rendezvous 
technology  has  also  evolved  with  small  spacecraft  development,  such  as  the  National 
Aeronautics  and  Space  Administration’s  (NASA)  DART  (NASA  2004),  Air  Force 
Research  Laboratory’s  (AFRL)  (Experimental  Satellite  System)  XSS-10  (Davis  2003) 
and  XSS-1 1  (AFRL  2005),  Naval  Research  Laboratory  (NRL)  SUMO  (Bosse  2004),  and 
Defense  Advanced  Research  Projects  Agency’s  (DARPA)  Orbital  Express  (OE) 
(Kennedy  2008).  In  particular,  the  AFRL  Space  Vehicle  Directorate  at  Kirtland  Air 
Force  Base  in  New  Mexico  developed  the  XSS-1 1  in  order  to  exhibit  the  ability  for  a 
small  satellite  to  autonomously  plan  and  rendezvous  with  a  passive  and  cooperative 
Resident  Space  Object  (RSO)  in  low  Earth  orbit  (LEO).  XSS-10  was  a  simple  proximity 
mission  around  the  upper  stage  that  it  was  boosted  into  orbit  on.  XSS-1 1  demonstrated 
the  ability  to  change  orbits  and  intercept  other  nonrelated  objects.  The  use  of  micro¬ 
satellites  to  inspect,  service,  repair,  and  refuel  larger  spacecraft  is  a  long-term  goal.  The 
closest  the  XSS-1 1  approached  and  maneuvered  around  another  object  in  space  was 
approximately  500  meters.  In  addition,  DARPA’s  OE  Advanced  Technology 
Demonstration  Program  validated  the  technology  and  techniques  for  on-orbit  refueling 
and  reconfiguration  of  two  satellites.  The  mission,  conducted  in  2007,  performed  several 
autonomous  rendezvous  and  capture  scenarios,  including  component  exchange  and 
propellant  transfer  events.  The  existence  of  these  programs  demonstrates  that  there  is  a 
need  for  a  robust  and  effective  autonomous  close  proximity  control  algorithm  for 
multiple  small  spacecraft.  Still,  even  with  the  investment  in  the  aforementioned 
missions,  spacecraft  proximity  operations  that  have  any  dynamic  tasking  attributes  are 
currently  executed  with  humans  in  the  loop,  relying  heavily  on  an  extensive  operations 
staff  to  plan  and  oversee  the  maneuvers,  and  performing  maneuvers  while  in 
communication  with  ground  stations.  This  approach  is  extremely  cumbersome  with 
respect  to  time,  manpower  and  overall  operational  footprint. 
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From  a  historical  point  of  view,  one  of  the  first  attempts  to  recapture  a  satellite  in 
orbit  was  in  1992  with  Intelsat6  from  the  Space  Shuttle  (Broad  1992).  In  this  case,  jet 
thrusters  on  both  the  shuttle  and  the  satellite  were  repeatedly  fired  during  the  three-day 
operation  to  bring  the  two  spacecraft  to  a  complex  orbital  rendezvous.  The  ground 
controllers  of  Intelsat  had  slowed  the  rotation  of  satellite  to  0.67  revolutions  per  minute 
so  it  would  be  easier  to  grasp.  Still,  astronauts  in  spacesuits  tried  and  failed  twice  capture 
the  satellite  from  orbit  and  move  it  into  the  payload  bay  of  the  space  shuttle  Endeavour. 
Success  was  eventually  realized  on  the  third  try. 

The  vast  majority  of  previous  research  pertaining  to  operations  that  focused  on 
rendezvous  with  an  uncontrolled  rotating  object  has  only  addressed  the  translational 
challenges  (Jacobsen,  Lee,  Zhu  and  Dubowsky  2002;  Matsumoto,  Dubowsky,  Jacobsen, 
and  Ohkami  2003a  and  2003b).  This  also  involved  segmented  maneuvers  and  only 
matching  the  translational  position  and  velocity  of  a  point  mass  with  that  of  a  docking 
point  on  a  tumbling  RSO.  Nolet,  Kong,  and  Miller  used  a  simple  glideslope  algorithm 
for  the  purpose  of  attempting  to  generate  trajectories  for  the  two-dimensional  (2D)  case 
with  hopes  of  implementation  on  the  Massachusetts  Institute  of  Technology  SPHERES 
testbed  (2004).  Huntington  (Huntington  and  Rao  2008b)  and  Singh  (Singh  and  Hadaegh 
2001)  independently  researched  proximity  operations  but  with  cooperative  agents. 
Henshaw  mentions  optimal  rendezvous  with  a  tumbling  object,  but  there  is  no  formal 
presentation  of  the  results  using  the  Minimum  Principle  (MP)  (Henshaw  2003). 
Henshaw’ s  method  is  a  variation  trajectory  planner  capable  of  planning  large  timescale 
maneuvers  (on  the  order  of  10,000+  seconds  to  several  years),  therefore  there  is  no 
discussion  of  the  computational  times  involved  to  obtain  a  solution. 

The  optimal  satellite  reorientation  problem  alone  is  of  general  interest  to  many  in 
the  field  of  aerospace  engineering.  Its  interest  increases  as  we  start  to  consider 
autonomous  rendezvous  and  close  approach  mission  scenarios  that  couple  both  attitude 
and  translational  dynamics  with  the  desire  to  match  both  translational  and  angular  motion 
in  preparation  for  docking.  The  available  literature  on  optimal  and  efficient  spacecraft 
reorientation  alone  is  extensive  (Junkins  and  Turner  1986;  Bilimoria  and  Wie  1993; 
Vadali  and  Junkins  1984).  Many  civilian  and  military  space  missions  need  to  have  agile 
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attitude  maneuver  capability.  For  instance,  Tactical  Satellite  (TacSat)  3  was  intended  to 
demonstrate  responsive  delivery  of  information  to  operational  users  (Davis  and  Straight 
2006).  Due  to  the  satellite’s  Low  Earth  Orbit,  the  timeline  for  tasking,  slewing  and 
disseminating  data  is  greatly  reduced.  Other  challenges  include  the  fact  that  the  tasking 
can  be  modified  at  any  moment  until  after  a  short  period  of  time  with  the  ground  station 
who  is  uploading  the  tasking  and  the  idea  that  TacSat  must  autonomously  slew  to  the 
target,  collect  and  process  the  data  and  then  down  link  the  data  directly  to  the  customer 
who  is  not  collocated  with  the  ground  station.  Current  real  time  feedback  controls  are  not 
optimized  for  minimum  time  (Vadali  and  Junkins  1984).  These  challenges  lend 
themselves  to  the  need  for  the  ability  to  rapidly  generate  feasible  trajectories  that  are 
optimized  for  minimum  time. 

The  preceding  examples  bring  about  a  very  interesting  concept:  what  if  only  a 
feasible  solution  is  desired,  with  little  or  no  regard  for  maneuver  time  (aside  from  a 
maximum  final  time)  or  power  used,  where  mission  success  is  the  only  criterion 
evaluated?  Is  the  new  metric  for  optimization  based  on  computational  time  or 
robustness?  Such  is  the  case  for  the  reorientation  of  a  Tactical  Satellite  (TacSat)  (Davis 
and  Straight  2006).  The  desire  is  to  perform  the  mission  within  certain  operating 
constraints.  If  the  satellite  happens  to  reorient  in  minimum  time,  nothing  will  be  gained 
(as  long  as  all  other  constraints  on  the  dynamics  are  met)  because  the  target  or  the 
downlink  may  not  be  in  view,  or  the  elevation  mask  to  the  target  may  not  be  cleared.  For 
certain  missions,  due  to  target  spacing,  access  times  and  general  Concept  of  Operations 
(CONOPS),  a  hurry-up  and  wait  mentality  will  not  benefit  the  mission  as  much  as  rapidly 
generating  the  trajectory  for  the  satellite  to  be  dynamically  tasked  at  the  latest  possible 
instant  (not  orbits  ahead  of  time).  In  order  for  the  satellite  to  perform  the  maneuver,  with 
potentially  target  locations  dynamically  changing  up  until  the  point  the  link  to  the 
groundstation  is  broken,  a  feasible  trajectory  needs  to  be  calculated  in  the  minimum 
amount  of  time  possible  before  the  first  target  acquisition.  In  essence,  there  is  no  time  to 
wait  for  a  detailed  optimal  solution  to  be  computed,  but  the  implementation  of  the  best 
optimized  solution  at  that  time  is  of  utmost  importance.  This  concept  will  be  revisited 
throughout  the  different  subsections  of  this  dissertation  where  the  rapid  generation  of  the 
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solution  and  robustness  of  always  having  a  feasible  trajectory  to  follow,  is  considered 
more  important  than  the  perfonnance  index  (PI),  such  as  minimum  time,  minimum 
energy,  etc.,  to  be  optimized. 


The  entire  mission  of  rendezvous  and  docking  can  be  broken  down  into  several 
components  and  subcomponents.  Figure  1,  modified  from  Fehse  (2003),  illustrates  the 
various  segments  of  the  mission. 


Docking 

Structural  Connection 


Capture 


_ t _ 

Close  Range  Rendezvous 

Final  Approach 
Approach  to  capture  point 
Achievement  of  final  conditions 

Closing 

Reduction  of  relative  distance  to  target 
Preparation  for  final  approach 


_ t _ 

Far  Range  Rendezvous 

Transfer  from  Phasing  into  close 
vicinity  of  target,  relative  navigation 

t 

Phasing 

T 

Launch 

Figure  1 .  Breakdown  of  rendezvous  and  docking  mission  elements  (After  Fehse  2003). 

The  work  in  this  dissertation  revolves  around  the  desire  to  efficiently  perform 
close  range  rendezvous  maneuvers  to  a  tumbling  spacecraft  for  the  purpose  of  close 
inspection  or  docking.  The  problem  of  docking,  contact  and  grappling  once  the  chaser 
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satellite  is  in  position  is  currently  being  explored  by  NRL  (Creamer  2007;  Tasker  and 
Henshaw  2008).  Among  others,  research  by  McCamish  (2007)  addresses  the  problem  of 
efficiently  and  safely  taking  a  satellite  from  far  range  into  close  proximity,  or  closing,  to 
an  RSO,  albeit  a  nontumbling  one.  Therefore,  this  research  fits  in  the  realm  of  close 
approach  trajectory  generation,  or  more  specifically  final  approach,  for  the  purpose  of 
close  inspection  or  docking  with  a  tumbling  RSO.  Launch  and  Phasing  (maneuvering  to 
match  orbital  elements  and  far  approach)  are  assumed  to  have  already  taken  place. 
Specifically,  adapting  the  tenninology  set  forth  by  Nolet  (2007),  this  is  subset  of  the 
tenninal  phase,  starting  ~10m  from  the  RSO  up  until  docking  contact  has  been  initiated. 
In  concert  with  previous  research  in  this  area,  it  is  not  the  intent  to  advance  navigation 
techniques  for  determining  the  attributes  of  a  tumbling  RSO,  therefore  it  assumes  that 
exact  knowledge  of  the  RSO  states  are  available.  Complementary  research  efforts  are 
looking  at  the  problem  of  detennining  the  states  of  an  RSO  (Jasiobedski,  Greenspan  and 
Roth  2001;  Tsuda  and  Nakasuka  2003).  An  exhaustive  list  of  guidance,  navigation  & 
control  (GNC)  architecture  efforts  related  to  rendezvous  can  be  obtained  in  Nolet  (2007). 
The  specific  research  addressed  here  currently  has  direct  application  to  a  specific  type  of 
problem  scenario  posed  by  Nolet  (2007)  to  rendezvous  with  a  target  that  is  “drifting  and 
tumbling”  as  defined  below: 

Target  drifting  or  tumbling:  it  is  able  to  communicate  its  states,  but  has  no 
control  authority  on  its  displacements  and  attitude  (e.g.,  docking  with  a 
fuel-depleted  target  with  no  reaction  wheels). 

Targets  that  are  fully  cooperative  or  able  to  communicate  and  control  attitude 
information  would  also  be  candidates  as  considered  by  Romano  (Romano  and  Hall  2006; 
Romano,  Friedman  and  Shay  2007).  This  research  would  be  coupled  with  advanced 
navigation  methods  to  attempt  to  explore  scenarios  where  there  is  no  information  sharing 
between  satellites  as  listed  below  (Nolet  2007): 

Target  not  cooperating: 

•  the  target  has  no  control  authority,  and  no  communication  is 
occurring  (e.g.,  docking  with  a  dead  satellite); 
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•  the  target  is  actively  trying  to  escape  the  chaser,  and  no 
communication  is  occurring  (e.g.,  military  applications). 

Regardless  of  the  specific  application  and  overall  GNC  architecture,  the  trajectory 
generation  methods  will  assume  to  have  state  information  regarding  the  attitude  and 
position  of  the  tumbling  RSO. 

This  dissertation  is  organized  as  follows:  Chapter  I  contains  a  summary  of  past 
and  current  research  on  this  topic  and  an  outline  of  several  concepts  employed  in  the 
dissertation.  Chapter  II  focuses  on  the  6  DoF  spacecraft  rendezvous  problem  while 
formulating  the  dynamic  equations  for  matching  points  of  interest  on  two  separate 
spacecraft,  one  of  which  is  controlled,  the  other  of  which  is  tumbling.  The  optimal 
control  problem  is  addressed,  resulting  in  deriving  the  adjoint  equations  and 
transversality  conditions  needed  to  verify  optimality.  Solutions  are  obtained  using  a 
direct  method  and  are  verified  using  indirect  shooting  methods  based  on  the  MP 
(Pontryagin,  Boltyanskii,  Gamkrelidze,  and  Mishchenko  1964).  The  purpose  of  this 
development  is  to  provide  a  baseline  solution  to  the  formulated  close  approach  problem 
that  has  not  yet  been  addressed  in  research.  Chapter  III  expands  on  the  results  of  Chapter 
II  increasing  the  complexity  of  the  problem.  The  same  approach  is  taken,  but  the 
translational  actuators  (thrusters)  are  body  mounted,  the  inertia  matrix  of  the  tumbling 
object  is  no  longer  taken  to  be  identity  and  the  end  constraints  require  not  only  matching 
position  and  velocity  of  the  spacecraft  docking  points,  but  also  their  attitude  and  angular 
velocity.  Chapter  IV  introduces  the  concept  of  Inverse  Dynamics  in  the  Virtual  Domain 
(IDVD)  in  more  detail  as  it  applies  to  spacecraft  reorientation,  a  key  enabler  for 
operations  with  respect  to  a  tumbling  object  that  has  not  been  addressed  so  far.  Chapter 
V  develops  a  method  for  rapid  trajectory  generation  based  on  IDVD  and  applies  to  the 
spacecraft  rendezvous  problem.  Chapter  VI  contains  analysis  and  simulations  of  the 
trajectories  generated  using  the  methods  developed.  This  includes  comparing  the 
performance  to  simplified,  2D  cases  that  exist  in  literature  and  evaluating  the  overall 
effectiveness  when  subjected  to  the  tracking  challenges  of  current  realistic  systems.  The 
chapter  concludes  with  ideas  for  implementation  in  closed-loop  architecture  that  would 
exploit  the  calculation  and  robustness  advantages  of  IDVD  methods. 
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A.  OPTIMAL  CONTROL  PROBLEM  FORMULATION  FOR  CLOSE 

RENDEZVOUS  AND  DOCKING 

Previously,  the  optimal  control  for  rendezvous  with  a  tumbling  object  has  only 
been  studied  in  2D  (Ma,  Ma  and  Shashikanth  2007)  with  very  little  analysis.  Studying 
the  optimal  solution  of  a  given  problem  also  entails  verifying  its  optimality.  This  is  done 
through  rigorous  analysis  of  the  states,  costates  and  how  they  relate  to  the  MP.  This 
research  formulates  the  three-dimensional  (3D)  problem  and  utilizes  a  direct  collocation 
method  (specifically  a  pseudospectral  method)  to  obtain  a  solution,  which  is  then 
validated  using  an  indirect  method  based  on  the  MP.  Although  direct  collocation 
methods  enable  the  generation  of  an  optimal  solution  (while  indirect  methods  have  a  very 
limited  radius  of  convergence),  they  are  extremely  computationally  intensive  and  require 
a  significant  amount  of  time  to  converge  (Yakimenko,  Xu  and  Basset  2008,  Boyarko, 
Yakimenko  and  Romano  2009a  and  2009b).  Therefore,  other  rapid-trajectory  generation 
techniques  are  developed,  while  using  the  optimal  solution  as  a  baseline  for  comparison. 

B.  INVERSE  DYNAMICS  IN  THE  VIRTUAL  DOMAIN 

1.  The  Concept  of  Inverse  Dynamics  in  the  Virtual  Domain 

On  top  of  having  an  extremely  large  computational  time,  direct  collocation 
methods  may  also  converge  to  sub-optimal  solution  if  the  number  of  nodes  is  not 
sufficient  (Yakimenko  et  al.  2008;  Boyarko  et  al.  2009a  and  2009b).  In  addition  to  that, 
the  optimal  solution  does  not  have  an  analytical  representation,  which  may  pose  problems 
when  trying  to  implement  it  using  a  feed-forward  scheme  of  suggested  control  commands 
(Yakimenko  et  al.  2008;  Boyarko  et  al.  2009a  and  2009b).  Furthermore,  the  placement  of 
nodes  by  a  direct  collocation  method  that  represent  the  solution  are  not  flexible  can  cause 
problems  when  trying  to  interpolate  complex  control  solutions  that  are  not  necessarily 
continuous  (Hurni  2009). 

This  research  pursues  another  approach  exploiting  the  general  idea  of  the  direct 
optimization  methods  of  calculus  of  variations  together  with  an  inverse  dynamics 
approach.  In  particular,  polynomials  are  used  as  basis  functions  to  generate  trajectories 
that  can  be  traversed  according  to  a  computed  analytic  speed  profile.  An  abstract 
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argument  is  used  that  allows  the  trajectory  to  be  decoupled  in  space  and  time,  while 
always  resulting  in  a  trajectory  that  will  reach  the  desired  endpoint  conditions.  The 
resulting  quasi-op timal  trajectory  solution  is  capable  of  being  generated  (and  updated) 
rapidly  because  of  the  reduction  in  the  number  of  varied  parameters  due  to  the  restriction 
on  the  trajectory  structure  by  specifying  a  polynomial  basis.  Although  this  method  lacks 
some  flexibility  due  to  the  predefined  structure,  it  provides  a  feasible  solution  that 
satisfied  the  endpoint  constraints  on  the  trajectory  at  every  iteration,  even  when  the  initial 
conditions  change  due  to  disturbances  or  delays.  Specific  applications  include  scenarios, 
where  derivative  conditions  on  beginning  and  ending  states  need  to  be  met,  such  as 
tracking  missions,  docking  missions  and  other  missions  that  are  not  necessarily  rest-to- 
rest  maneuvers.  It  also  has  application  in  areas,  where  computational  time  and  the  ability 
to  immediately  employ  the  best  available  solution  at  any  given  time,  is  heavily  weighted 
as  part  of  the  overall  PI  of  the  solution. 

2.  Attitude  Trajectory  Generation 

The  goal  of  this  research  is  to  provide  a  method  to  determine  a  feasible  attitude 
trajectory  solution  that  meets  endpoint  requirements  and  dynamic  constraints  while 
performing  a  good  overall  maneuver  relatively  to  a  given  cost  function.  Also,  the  method 
should  work  for  any  terminal  conditions  including  nonrest  to  nonrest  maneuvers.  The 
major  requirement  is  that  the  method  must  focus  on  providing  a  rapid,  potentially  real 
time  solution  as  opposed  to  off-line  computations  even  if  it  requires  some  sacrifices  in 
optimality.  This  methods  developed  here  are  coupled  with  methods  developed  in  a  later 
section  for  the  translational  trajectory  generation,  providing  the  complete  solution  needed 
for  close  approach  and  docking  with  a  tumbling  RSO. 

3.  Translational  Trajectory  Generation 

Trajectory  generation  by  IDVD  has  been  used  in  the  past  for  aeronautical 
applications  where  the  only  two  frames  that  needed  consideration  were  the  inertial  and 
body  frames  (Yakimenko  et  al.  2008;  Yakimenko  and  Siegers  2009).  This  research 
utilizes  the  concept  of  IDVD  in  the  orbital  frame  fixed  at  the  RSO  (an  intermediate  frame 
between  the  inertial  and  body  frame)  and  exploiting  the  Clossey-Wilshire  Equations  of 
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relative  motion.  This  allows  constructing  the  problem  to  analytically  set  the  position  and 
rate  values  of  the  states  in  the  specific  frame  of  interest,  ensuring  a  maneuver  that  will 
perform  the  desired  close  approach  maneuver  with  respect  to  the  RSO,  allowing  for  rapid 
trajectory  reshaping  or  reoptimizing  and  reducing  overall  relative  error.  The  use  of  a 
speed  factor  of  a  specific  form  (having  an  analytical  expression  for  the  integral  of  the 
inverse  function),  also  allows  flexibility  in  setting  the  node  points  at  specific  time 
intervals.  This  new  method  of  rapid-trajectory  generation  allows  the  computation  of  an 
optimized,  feasible,  safe  trajectory  for  spacecraft  rendezvous  requiring  significantly  less 
computational  time  than  current  methods. 

C.  SCOPE  OF  THE  DISSERTATION 

This  dissertation  advances  the  body  of  knowledge  with  respect  to  guidance  of  a 
spacecraft  to  rendezvous  and  dock  with  a  tumbling  object.  In  particular,  the  main  new 
contributions  of  this  dissertation  are: 

1.  The  analytical  fonnulation  of  the  optimal  control  problem  is  presented  with 
respect  to  a  spacecraft  perfonning  a  controlled  approach  in  3D  to  a  tumbling 
object  for  the  purpose  of  rendezvous  and  docking.  A  baseline  optimal 
solution  is  presented  using  pseudospectral  methods  and  verified  the  necessary 
conditions  for  optimality  based  on  the  Minimum  Principle  coupled  with  an 
indirect  shooting  method. 

2.  An  analytic  quaternion  trajectory  is  formulated  using  a  5th  and  7th  order 
Bezier  Polynomial.  Derivatives  are  formulated  and  applied  in  the  virtual 
domain  providing  an  attitude  trajectory  generation  technique  that  can 
automatically  match  specified  derivative  values  at  the  endpoints  while  varying 
the  speed  across  the  spatial  trajectory.  This  fonnulation  is  then  applied  for 
rapid  trajectory  generation,  significantly  decreasing  the  computational  time. 

3.  The  analytic  quaternion  trajectory  approach  based  on  Bezier  polynomials  is 
coupled  with  a  translational  polynomial  scheme  to  create  a  novel  trajectory 
generation  scheme  for  the  close  approach  with  a  tumbling  object.  The  method 
is  extremely  flexible  allowing  for  real  time  reshaping  of  the  trajectory  based 
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on  dynamically  changing  end  constraints.  Also,  by  incorporating  a  speed 
factor  of  a  particular  form,  this  method  can  provide  continuous  mapping  from 
the  virtual  and  time  domains.  Analysis  shows  that  a  trajectory  can  be 
calculated  that  is  within  -10%  of  the  optimal  trajectory,  but  only  requiring  a 
fraction  (-0.5%)  of  the  computational  time. 

4.  The  IDVD  method  is  incorporated  in  different  GNC  architecture  schemes  to 
demonstrate  its  flexibility.  An  IDVD  reshaping  method  is  created  and 
implemented  in  a  real  time  rapid- trajectory  generator  that  incorporates  the 
most  current  knowledge  of  the  RSO  states.  A  hybrid  method,  using  a 
combination  of  recalculating  and  reoptimizing  the  solution  as  well  as 
reshaping  based  on  current  RSO  state  knowledge,  is  developed  and 
demonstrated. 
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II.  FORMULATION  AND  ANALYSIS  OF  MATCHING  POINTS  OF 
INTEREST  IN  TWO-SPACECRAFT  FOR  OPTIMAL 

RENDEZVOUS 

From  the  theoretical  standpoint,  this  chapter  elaborates  on  the  previous  work  by 
Ma  (Ma  et  al.  2007),  who  has  studied  the  minimum-control  effort  problem  for  a  planar 
rendezvous  to  a  tumbling  object  (with  only  three  states,  coordinates  x  and  y,  and  heading 
angle  6),  neglecting  any  path  constraints  and  relative  motion  dynamics  pertinent  to 
proximity  space  operations.  The  Sakawa-Shindo  algorithm  was  used  for  calculating  the 
optimal  control.  This  chapter  expands  the  scope  by  taking  into  account  the  proximity 
motion  dynamics,  considering  the  full  six  DoF  model,  and  detennining  both  the 
minimum  time  and  the  quadratic  fonnulation  of  minimum-control  (energy)  effort  solution 
of  the  rendezvous  of  two  satellites.  It  also  features  a  comparison  of  the  solutions  obtained 
using  one  of  the  prominent  direct  collocation  methods  with  the  truly  optimal  solutions 
obtained  using  the  MP.  Another  section  of  this  dissertation  takes  an  even  further  step  and 
considers  three  different  Pis,  adding  additional  constraints  to  match  tenninal  attitude  and 
angular  dynamics,  along  with  position  and  velocity. 

A.  TWO-SPACECRAFT  RENDEZVOUS  MODELING  AND  OPTIMIZATION 

PROBLEM  FORMULATION 

This  section  develops  a  model  of  target-chaser  rendezvous.  Figures  2  and  3  show 
a  graphical  representation  of  the  problem.  The  center  of  the  orbit  frame  is  fixed  to  the 
center  of  mass  (CM)  of  the  tumbling  RSO.  The  x  axis  points  toward  the  Zenith.  The  y 
axis  lies  along  the  velocity  vector  of  the  RSO  (assuming  circular  orbit),  and  the  z  axis  lies 
along  the  orbit  normal  of  the  RSO.  We  start  from  the  arbitrary  relative  position  (Figure  2) 
and  attempt  to  bring  two  spacecraft  together  for  docking  (Figure  3). 
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Figure  2.  Depiction  of  the  spacecraft  rendezvous  problem. 


Figure  3.  The  desired  final  state  of  two  spacecraft  system. 
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The  set  of  notations  to  be  used  in  this  model  (some  of  which  are  shown  on  Figures 
2  and  3)  include: 


m 

df ,  i  =  x,y,z 


df ,  i  =  x,y,z 


x,y,z 


x,y,z 


cof ,  i  =  x,y,z 


CO? ,  i  =  x,y,z 


o  c 

cot  ,  i  —  x9y9z 


O  T 

cot  ,  i=x9y,z 


q? ,  i  =  1,2, 3, 4 


the  mass  of  the  chaser  spacecraft; 

the  position  vector  of  the  target  docking  point  with  respect  to  the 
target  CM  expressed  in  the  target  body  frame; 

the  position  vector  of  the  chaser  docking  point  with  respect  to  the 
chaser  CM  expressed  in  the  chaser  body  frame; 

the  x,  y,  and  z  position  of  a  chaser  spacecraft  CM  in  the 
Hill’s  coordinate  frame; 

the  x,  y,  and  z  component  of  the  velocity  of  chaser  spacecraft 
center  of  mass  in  the  Hill’s  coordinate  frame; 

the  x,  y,  and  z  component  of  the  angular  velocity  of  the  chaser 

spacecraft  with  respect  to  the  inertial  frame,  expressed  in  the 
chaser  spacecraft  principal  body  coordinate  frame; 

the  x,  y,  and  z  component  of  the  angular  velocity  of  the  target 
spacecraft  with  respect  to  the  inertial  frame,  expressed  in  the  target 
spacecraft  principal  body  coordinate  frame; 

the  x,  y,  and  z  component  of  the  angular  velocity  of  the  chaser 

spacecraft  with  respect  to  the  orbital  frame,  expressed  in  the  chaser 
spacecraft  principal  body  coordinate  frame; 

the  x,  y,  and  z  component  of  the  angular  velocity  of  the  target 

spacecraft  with  respect  to  the  orbital  frame,  expressed  in  the  target 
spacecraft  principal  body  coordinate  frame; 

the  components  of  the  quaternion,  representing  rotation  from  the 
orbital  to  chaser  body  frame; 
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c/l ,  i  =  1,2, 3, 4  the  components  of  the  quaternion,  representing  rotation  from  the 

orbital  to  target  body  frame; 

Q  the  angular  rate  of  the  orbital  frame  with  respect  to  the  inertial 

frame. 

Using  these  notations,  the  dynamics  of  the  two  systems  can  now  be  described  as 
follows  (Figure  2).  The  translational  kinematics  and  dynamics  of  a  chaser  spacecraft  in 
the  orbit  frame  centered  at  the  target  vehicle  are  presented  by  Hill’s  equations  (Vallado 
and  McClain  2007): 


x  =  — (2Qv  +  3 Cl2x  +  fx) 
m 

y  =— (-2Qx + f  )  (1) 

m 

z  =  -{-Orz  +  f„) 
m 


where  fx  ,  f  and  f_  are  the  applied  forces  (controls)  expressed  in  Hill’s  frame. 

The  rotational  dynamics  of  the  chaser,  described  by  the  vector  equation 
(Greenwood  1987;  Wie  1988) 

Ictdc+(ocxIctoc  =T  (2) 


expands  into  the  scalar  quantities: 


ox 


(/£-4)©>f+7? 


CD,.  =- 


I 


22 


of 


33 


(3) 
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In  Equations  (2),  (3)  I  =  diag([I^ ,  /£ ,  /53 ])  is  the  inertia  matrix  along  the 
principal  axes,  (oc  =  [mcx  ,ry^,ruf]T ,  and  T  =  [T  ,  7)  ,  7)]  '  is  a  vector  torques  (symbol  T 
denotes  transposition). 

Similarly,  for  the  target  (we  have  no  control  of)  the  rotational  dynamics  is  given 
by  the  vector  equation 

iW+G^xlW  =0  (4) 

or  its  scalar  form 


Defining  °R  as  the  rotation  matrix  to  convert  from  the  body  frame  of  the  target 

{T}  to  the  orbital  frame  {0}  and  °c  R  to  convert  from  chaser  body  frame  {C}  to  the 
orbital  frame  {O},  we  can  define  the  angular  velocity  of  each  object  ( a  =  {T,C })  with 
respect  to  the  orbital  frame  expressed  in  the  orbital  frame: 

r°_ 

R  <  -  o  (6) 

o.fz  |_Q 


Rotation  matrices  „R  are  constructed  using  components  of  the  corresponding 
quaternion  as  follows  (Greenwood  1987;  Wie  1988): 

+q‘l-q?-q?  2(q°q*,  -q°q°t)  +  9?94") 

°R.  2  (q-q°  +  q-q;)  qf-qf+qf-qf  2  (q°2q‘  -  q?q", )  (7) 

2(q?q“~qW)  2  (qiql+q°q“)  -  qf  ~  qf  +?f  . 
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Quaternions  are  a  parameterization  for  spacecraft  attitude  transfonnations  based 
on  the  fact  that  the  general  displacement  of  a  rigid  body  with  one  fixed  point  is  a  rotation 
about  a  fixed  axis:  the  eigenaxis.  The  definition  of  the  quaternion  is: 

px  sin( cr  /  2) 

A  sinter  /  2) 

q=  (8) 

A  sin(cr/2) 
cos(cr/2) 


where  cr  is  the  scalar  value  of  the  rotational  displacement  about  the  eigenaxis  and  p  is 
unit  vector  that  describes  the  direction  of  the  eigenaxis.  The  parameters  of  quaternions 
are  propagated  according  to  the  known  relation  (Greenwood  1987;  Wie  1988): 


ql 

i 

ql 

~2 

a:\ 

O  ..a 

—  CO, 
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COy 

-V“ 


O  a 

CO, 

1 

o 

v8 

^  £ 

o 

§ 

*  S5 

_ 1 
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0)x 

C>  a 
°>y 

A 

1 

O 

a 

*  £ 

0 

O 

CD, 

A 

O  a 
~  % 

0  a 

-  CD, 

0 

1 - 

£ 

(9) 


Equations  (1),  (3),  (5)  and  (9)  define  a  20-state  system  of  differential  equations  governing 
rendezvous  dynamics.  Combined  into  the  state  vector  x  these  states  are: 


x  [x5  y,  z,  x,  y,  z,  cox ,  co  ,  co,  ,  cox ,  co  ,  co,  ,qx  ,  A  -  A  -  A  -  q\  -  A  •  7'.  ■  7  i  I 


(10) 


The  governing  dynamics  assume  six  normalized  controls  that  can  be  used  by  the 
chaser  to  achieve  the  rendezvous  conditions: 


u 


/.  f,  f. 


f  f  f 

J  imax  J  vmax  J  z i 


T  T 

y  max  z  r 


(11) 


For  simplicity,  we  will  further  assume  /)max  =  \m  /  s2 ,  T.mm  =  \Nm ,  i  =  x,y,z. 
These  controls  are  three  normalized  components  of  a  translational  force  acting  on  the 
chaser  ft,  i  =  x,y,z,  expressed  in  the  Hill’s  coordinate  frame,  and  three  normalized 
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components  of  a  torque  allowing  to  change  chaser’s  attitude,  T.,  i  =  x,y,z  ,  expressed  in 
the  chaser’s  body  frame.)  All  six  controls  are  bounded:  -1  <  u  <  1 . 


Using  these  controls,  we  would  like  to  bring  the  two  spacecraft  from  some  initial 
conditions,  given  by  20  initial  values  of  states  xj(t0) ,  i  =  1, ...,  20 ,  to  docking-enabling 
conditions  described  by  matching  chaser’s  and  target’s  docking  station  final  positions  and 
velocity  vectors  as  described  by  the  following  equations.  These  conditions  for  positions 
and  velocities,  expressed  in  the  matrix  form  are  as  follows: 


„R 


*  M 
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1 

^  * 

1 _ 

X 
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- 
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= 

e2 
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_e3_ 

=  0 


and 


1 

*  ^ 
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^  * 

^3 
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A6_ 

(12) 


While  transitioning  to  docking-enabling  conditions,  we  would  like  to  minimize  and 
compare  two  different  performance  indices: 

‘f 

J  =  \fodt  (13) 

*0 


with 


/o=l 


(14) 


for  minimum  time  and 


f0  —  —  [ll^  +  +  U-^  +  (15) 

for  minimum  quadratic-control  (or  more  commonly  referred  to  as  minimum  energy) 
expenditure  (Henshaw  2003).  An  equivalent  formulation  to  the  minimum  time  cost  can 

also  be  expressed  as  J  =  tf,  which  employs  an  endpoint  cost  (Mayer)  as  opposed  to  a 
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running  cost  (Lagrange).  A  brief  statement  about  the  cost  function  (or  PI)  is  in  order. 
The  specific  cost  function  for  minimum  energy  was  chosen  to  include  all  actuators.  The 
fact  is,  it  can  contain  any  subset  of  the  controls  or  states.  The  choice  for  this  research  was 
to  keep  the  cost  function  as  generic  as  possible  (including  all  actuator  in  the  cost),  but  one 
can  easily  think  of  situations,  where  the  cost  function  might  be  truncated  to  only  include 
translational  actuators  (the  case  of  using  Control  moment  Gyroscopes  (CMG)  or  Reaction 
Wheels  (RW)  instead  of  thrusters  for  attitude  control).  The  specific  cost  function 
implemented  would  depend  on  the  details  of  the  mission,  system  and  the  discretion  of  the 
user,  regardless  of  the  method  applied  to  obtain  a  solution. 

B.  SYNTHESIS  OF  THE  OPTIMAL  CONTROL  USING  MINIMUM 

PRINCIPLE 

1.  Formulation  of  the  Optimal  Control  Problem 

We  start  from  the  general  formulation  for  the  Hamiltonian  of  the  system  with  the 
state  vector  x,  Equation  (10),  controls  vector  u,  Equation  (11),  and  running  cost  of 
Equations  (13)-(15): 

H  (k,  x,  u)  :=  /0  +  (k,  x)  (16) 

where  operator  (...)  on  the  right-hand  side  denotes  a  scalar  product  of  two  vectors,  and 
k  e  is  a  costate  vector  which  differential  equations  are  to  be  defined  later  in  this 
section. 

For  the  specific  system  of  Equations  (1),  (3),  (5)  and  (9),  with  the  running  cost 
from  Equation  (14)  (minimum  time  problem)  the  Hamiltonian  can  be  written  with  respect 
to  the  state  vector  x  defined  in  Equation  (10)  as: 
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H(k,x,  u)  :=  1  +  ^X4  +  A2x5  +AjX6  +— (A4(2Qx5  +3Q2Xj  +  u1)  +  A5(-  2Qx4  +  m2)  +  /16(-Q2x3  +u3 )) 

m v  ' 


+ 


( Ic -Ic\ 

1/22  133) 


XsXg  +  U4 


Ag  + 


(4-4) 


x7x9  +  u5 


^8  + 


\/ll  122 ) 


x8x7  + 1/6 


Ay 


22 


33 


(lT -IT)x  X  (lT -IT)x  X  ( f -f  ) 

\*22  ±33)a‘Wx\2  .  \  33  -'ll/A10-A12  .  (-'ll  122) 


4  + 


\  1  +’ 


IXnXio 


Aq  2 


22 


'33 


+4^((* 9  -(4  -4  -4  +  45)Q)*14  "(*8  -2(Vl5  -*13*16.)Q)*15  +  (X7  ~2(X13X15  +  *14*16  )  Q)  *16  ) 

+4  ^("(*9  "(4  "4  -4  +  4)Q)*13  +(X8  -2(Vi5  "  *13*16  )! Q)  *16  +  (*7  ~2(X,3X15  +  *14*16  )  Q)  *15  ) 
+4  \(~(X1  _2(x13x15  +x14x16)n)x14  +(x8  -2(x14x15  -x13x16)Q)x13  +(x9  -(4  “4  -4+4)Q)xi6) 
+4^(-(x7  -2(x13x15  +x14x16)Q)x13  -(x8  -2(x14x15  -x13x16)Q)x14  -(x9  -(4  “4  "4  +  4)q)*i5) 
"*"4  “1(^12  —  (*20  ~~  *n  *i8  “*"*19) £2) *is  —  (xn  —  2(x18x19  —  x17x20 )  Q)  x19  +  (*io  —  2(x17x19  +x18x20)Q^x20  j 
~4 8  (*12  —  (*20  ~ Xn  *is  *19 ) ^2^ *17  “*"(-^11  — 2(x18x19  —  x17x20)Q)x20  +  (*io  ~ 2(x17x19  +x18x20)Q)x19 

"*"4 “(_(*io  —  2(x17x19  +x18x20)Q^x18  +  (*n  —  2(x18x19  —  x17x20)Q)x17  +^x12  —  (-*20  —  x\i  —  *is  +  *19)  £2)  *20 
”*”2-20  ~z{^~{x\o  — 2(x17x19  +  x18x20 )Q)x17  —  (*11  —  2(x18x19  —  x17x20)Q)x18  —  ^x12  — (x20  —xxi  —  x18  +x19  jojx19  j 
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(the  possibility  of  a  singular  control,  when  X i+3(t)  =  0  ,  is  considered  in  Chapter  II. B. 3). 

Likewise,  developing  the  Hamiltonian  for  the  minimum  quadratic-control  cost 
function,  based  on  Equation  (15)  the  part  of  the  Hamiltonian  that  depends  on  the  controls 
is: 


H  (k, x, u)  i—  —  +  w-,2  +  u3  +  w4  +  w5  +  w6  j 

+-(X/lul+A5u2+X6u3)  +  -^A7u4+-^Asu5+--^X9u6 


Since  the  controls  enter  Hamiltonian  nonlinearly,  the  optimal  control  is  not  the 
bang-bang  control  anymore.  To  be  more  specific,  the  resulting  optimal  control  that 
minimizes  the  Hamiltonian  for  the  minimum  quadratic-control  cost  function  Equation 
(15)  (minimum  energy)  is  as  follows: 
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/ 


itl<_l 


tot 


c  ’ 


x. , 

-1  <  — 7X  <  1  for  /  =  4,5, 6,  where  k  =  i-  3 


X 


?+3 


‘tt 


>1 


Now,  the  differential  equations  for  costates  will  be  the  same  for  both  optimization 
problems  and  are  obtained  via: 


i  = 


'dH 
v  3x 


(22) 


The  first  six  adjoint  equations  are  given  by: 

X  =  -3X4Q2/«  1 ,  X,  =  0 ,  X  =  X6Q2m_1 , 

X4  =  — X  +  2Xflm  { ,  X  =  — X  -  2X4Qm_1 ,  X6  =  ~X  (23) 

The  next  six  adjoint  equations,  corresponding  to  the  states  7  through  12,  take  the  form  of: 


X  =  In  /33  V?  +  /z2  C/n  X*8  ~X3 

1 22  ^ 


33 


X4  ^  +  X5  2  ^-Y14  +  X6  2  ^"13 


(24) 


The  remaining  eight  adjoint  equations,  corresponding  to  the  states  13  through  20,  are  of 
the  form: 
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-  2x13  Q  -  2x16“Q  -  2x15~Q 


i0=-*0fafl*u-2xjixa-2x,ftcu) 

— /i]4  —  (  — —  (a'l6  —  X]3  X,  |  +  X|5)flJ 

- 4  ^-(2*15Qx14  +  (*s  -2(^15  -x13x16)Q)  +  2x16Qx13  +  2x13Q*16) 
-^6^(2^15^13  ~(X1  -2ixnxis  +xuxl6)D.)-2xl6DxH-2xnDx15) 


2|4  ~~  A3  ~  (2xi4  Q  +  (x9  ( 


_|  —  v2  _V2  +  v2 

\  a16  a13  a14  -r  a15 


^  O j  +  2x^  +  2x^  O 


1 


-314-(2x„nxM-2x16ax15-2x1!ar1() 


■As  ^(“(*7 -2(ws  +  *14*16 )Q)  +  2x16Qxi4  -2x15Q*13  +2x14Q*16) 

-^6^(2xi6a*i3  -(*8  -2(*14*15  -*13*16)Q)  +  2*15a*14  -2*14a*15) 


As  —  A3  2  (  2*i5Q*14  (*8  2(*14*15  *13*j6)Q^  +  2*15Q*14  2x13Q*16  ) 

-^4^(2*15^2*13  -2*16a*]4  +(*7  -2(*13*15  +x14x16)Q)-2x13Qx15) 


1 


-215-(2x13Qx,4-2x14nx,J  +  2xlsnx,6) 


■\6  —  |+2*132Q  +  2x142Q-(x9  -( 


1  1  2  2  2,2 

—  |x0  -|*16  -*13  -*14  +*15 


)q)  +  2x152Q 


Ae  =_A3^ (-2*16^14  -2*13f2*i5  +(*7  -2(*13*15  +  xHxl6)n)-2xl6Qx.l4) 

-^4i(2*16a*13+2*16a*13  +(*8  -2(x14x15  -*13*16)Q)-2*14a*15) 

-2j5  —  |2*132Q  +  2*142Q-(*9  -(*f6  -*j23  -x,24  +*j25)Qj  +  2*162Qj 


3ll2(2x„fix-14  2x13nx14+2x,!£2xls) 


(25) 


(The  adjoint  equations  for  A17,...,A20  can  be  derived  by  swapping  the  associated  chaser 
states  and  costates  with  the  RSO  variables  (ex.  A3  — »  A7>*i3  xn>  A3  A7>etc-))- 
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For  this  problem  formulation,  it  should  also  be  noted  that  the  third  and  the  sixth 
costate  equations  are  decoupled.  This  feature  will  be  addressed  further.  The  next  step  is 
to  examine  the  transversality  conditions,  which  with  account  to  the  following  coupled 
tenninal  variations,  followed  directly  from  Equation  (12).  The  six  scalar  equations  for 
the  components  of  the  vector  e  =  [e1,...,e6]'  represent  the  matching  of  position  and 
velocity  of  the  docking  points  in  the  three  components  of  the  orbital  coordinate  frame  at 
the  time  t  f .  If  positions  of  the  docking  points  are  matched  in  the  orbital  frame,  then 

matching  velocity  in  the  orbital  frame  will  lead  to  matching  velocity  in  the  inertial  frame 
as  well,  since  the  transport  theorem  will  have  the  same  effect  for  coincident  points. 

Let  variable  E  represent  an  endpoint  cost  that  is  set  to  zero  for  both  the  minimum 
time  and  minimum-control  scenarios.  Introducing  an  auxiliary  function,  E ,  several  more 
necessary  conditions  for  optimal  control  (i.e.,  transversality  conditions)  can  be  taken  into 
account  as  (Yan,  Fahroo  and  Ross  2002): 

E  =  E  +  j^ua  (26) 

(=1 


f 

v 


—  \T 

SET 

dx 


(27) 


The  first  six  equations  result  in  expressions  that  only  contain  the  parameters  u/ , 

i  =  1 . 6  : 

Tj  ( tj )  =  —ol ,  (tf )  =  —o2 ,  T,  (tj- )  =  — 1>3 , 

4,  (tf  )  =  -°4  >  A  (tf  )  =  -o5  ,\(tf)  =  -o6  (28) 

These  six  equations  are  then  properly  substituted  into  the  remaining  equations  to  give  the 
14  equations  for  the  end  conditions,  described  as  =  4]',  that  must  be  satisfied, 

along  with  the  six  docking  states  for  position  and  velocity  given  by  Equation  (12)  at  t  , , 
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and  the  twenty  known  initial  conditions.  For  example,  for  the  value  at  tf  of  the  7th  costate 
in  relation  to  the  first  transversality  condition  is: 


4  (^/)  (20l3*14  +*15*16)*23  2(xi3X[5  X14X16^X22 ) 


UA 


+ 

+ 


^2(xl3Xi5  Xl4X16)x2l  (xl6+xl3  X14  -^15 )-^23 ) ^ 

((•*16  "*"-*13  ^*14  ~X\5)X22  —  2C*13-*14  +  x\5X\s)X2\) 


(29) 


Then  substituting  from  Equation  (28)  leads  to: 

(Tl  (f/0  —  ^(ty)  — (2(x13X14  +  •*15'*16),*23  —  2(Xi3Xi5  —  X]4X16  )x22  )  (— ^4  (ty  )) 

+  (2(x13x15  - x14x16)x21  -  (4  +  x,2,  - xf4  - xf5 )x23 ) (~AS (t tf  ))  .  (30) 

+  ((4+4  -x,2,  -x2  )x22  -2(xl3x,4  +xl5Xl6)x2l)(— 44))  =  0 


The  same  is  done  for  the  remaining  13  equations. 

Therefore,  everything  is  in  place  to  solve  the  minimum  time  Minimum  Principle 
problem  converted  to  the  nonlinear  programming,  two  point  boundary  value  problem 
(TPBVP)  numerically.  Specifically,  given  the  dynamics  of  the  system  described  by 
Equations  (1),  (3),  (5),  (9),  and  adjoint  system  Equations  (23)-(25),  with  the  optimal 
controls  synthesized  as  Equation  (19)  or  Equation  (21),  with  the  initial  states  x;(t0) , 

i  =  l,. ..,20,  we  guess  on  the  values  of  the  initial  costates  44),  f  =  1,..., 20  and  the 
maneuver  time  t  ,  (2 1  variable  parameters),  in  attempt  to  zero  six  final  discrepancies  in 

matching  final  position  and  inertial  velocity  of  chaser’s  and  target’s  docking  stations 
(Equation  (12)),  together  with  fourteen  relations  resulting  from  the  transversality 
conditions  (Equation(26))  and  assure  that  (Y an  et  al.  2002) 

H(tf)  =  0.  (31) 


Therefore,  we  have  2 1  varied  parameters  and  2 1  conditions  to  satisfy. 


2.  Accounting  for  a  Possible  Path  Constraint 

For  the  close  rendezvous  problem,  we  should  obviously  take  into  account  one 

more  constraint.  This  constraint  is  that  the  CM  of  the  chaser  spacecraft  must  remain  at  a 
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distance  larger  than  some  minimum  distance  (a  “keep  out”  sphere  with  a  radius  r)  from 
the  CM  of  the  RSO,  which  is  coincident  with  the  origin  of  the  orbit  frame.  This  assures 
that  the  chaser  vehicle  will  not  pass  through  the  target  vehicle  in  order  to  reach  the 
docking  position.  Mathematically,  this  can  be  defined  as: 

h  =  (xf  +%2  +x] )-r2  >  0  (32) 

where  jc.  ,  /'  =  1,2,3  are  the  first  three  elements  of  the  state  vector  (9).  Furthermore,  while 
a  trajectory  is  on  a  path  constraint  h  =  0  ,  the  tangential  condition  must  also  be  satisfied 
(Pontryagin  et  al.  1964): 

h  =  —  =  xxx4  +  x2x5  +  x3x6  =  0  (33) 

dt 

Consequently,  the  Hamiltonian  should  be  augmented  by  another  term: 

H  (k,  x,  u)  :=  f0  +  (k,x)  +  juh  (34) 


where  ju  is  a  constant  and 


h  = 


d2h 
dt 2 


(35) 


(since  the  path  constraint  has  to  be  differentiated  with  respect  to  time  twice  before  the 
control  variables  appear  in  the  expression).  The  value  of  //  is  dictated  as  follows: 


/u  =  0  if  h>  0  =>  off  the  constraint  boundary 

ju>  0  if  /?  =  /?  =  /?  =  0  =>  on  the  constraint  boundary 

Therefore,  upon  first  contact  with  the  path  constraint,  the  costate  values  and 
Hamiltonian  will  be  discontinuous  (Bryson  and  Ho  1975). 

3.  Possibility  of  a  Singular  Control  for  a  Minimum  Time  Problem 

Upon  closer  inspection,  we  find  that  in  the  z-direction  translational  control  u2  is 
decoupled  from  all  other  controls.  From  Equation  (1)  we  have 
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x3  =  x6,  x6  =  m  '(-Q2.^  +  w3) 


(37) 


or 


and  Equation  (23)  yield 


x3  +  m  Li  x3  =  m  u3 


(38) 


Aj  =  m  1Q.2  A6,  A6  =  —Ai 


(39) 


or 

X6  +  m~1Q2A6  =  0.  (40) 

Now,  for  the  minimum  time  problem,  the  optimal  control  is  defined  as  in  Equation  (19), 
hence  for  u3  we  will  have 


f  1,  4<o 

u3=< 

3  l-i,  4>o 

With  account  of  Equation  (40),  we  can  state  that 

w3  (f)  =  f(Ai(t0),Ab(toj). 


(41) 


(42) 


Defined  by  the  natural  frequency  of  Equation  (40)  the  control  u3  can  switch  from 

^3 max  =  1  to  u3mm  =  —1  or  vice  versa  every  7TyfmCl  1  seconds.  Moreover,  the  solution 

to  Equation  (40)  with  the  initial  conditions  defined  by  \(t{))  and  X6 (tQ )  =  —Aj (tQ )  can 
be  found  analytically: 

A 6(t)  =  -A3(t0)m°'5Q.~1  sm(m  0 5Cit)  +  A6(t0) cos(m  0 5Cit) .  (43) 

The  roots  of  this  equation  are  defined  as 
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t  =  tan 


i 


m°- 5Q -1 . 


(44) 


Mo}m-°-5  n 
MO  / 


From  the  other  hand,  the  extremals  of  Equation  (43)  are  achieved  at 


t*  =  -  tan  1 


v  \  (to  ) 


0.5/^,- 

m  £  2 


J 


(45) 


As  seen  from  Equations  (43)-(45),  the  only  possibility  for  the  singular  control  would  be 
when 


A3(t0)  =  0  and  A 6(t0)  =  0 .  (46) 

In  this  case,  A 6(t)  =  0  and  the  optimal  control  cannot  be  defined  from  Equation  (41),  but 
requires  more  rigorous  analysis  (Bryson  and  Flo  1975). 

C.  METHODOLOGY  FOR  OBTAINING  A  SOLUTION  AND  CHECKING 

THE  OPTIMALITY 

This  section  presents  the  methodology  for  obtaining  and  verifying  optimal 
solutions  for  two  problems  posted  in  Chapter  II.  A.  Despite  the  fact  that  the  structure  of  an 
optimal  control  was  defined  analytically  (in  the  previous  section),  it  is  clear  that  there  is 
no  way  one  could  solve  this  problem  using  direct  shooting  approach  having  some 
arbitrary  values  of  varied  parameters.  No  matter  what  optimization  algorithm  to  solve  the 
TPBVP  would  be  used,  the  numerical  solution  will  diverge.  That  is  where  direct  methods 
of  calculus  of  variations  become  useful.  In  what  follows,  this  section  introduces  a  specific 
rendezvous  scenario  with  the  concrete  specific  numerical  values  to  play  with.  Next,  it 
describes  a  procedure  of  obtaining  quasi-optimal  numerical  solutions  for  each  of  two 
optimization  problems  using  direct  collocation  methods,  also  known  as  pseudospectral 
methods.  Finally,  a  methodology  of  using  this  solution,  which  is  very  close  to  the  true 
optimal  one,  to  address  the  problem  using  direct  shooting  method  for  the  TPBVP 
formulated  in  Chapter  II.B.l  is  introduced. 


27 


1.  Defining  a  Rendezvous  Scenario 

This  specific  scenario  to  play  with  initializes  the  chaser  CM  starting  at  a  distance 
of  5  meters  from  the  target  CM  with  the  target  having  an  initial  angular  velocity  of  0.25 
rad/s  in  both  y  and  z  body  axes  without  losing  generality.  The  body  coordinate  frames  of 
each  spacecraft  and  the  orbit  frame  are  assumed  to  be  coincident  with  the  inertial  frame  at 
the  beginning  of  the  simulation.  The  chaser  docking  point  is  located  at  [-0.25,  0,  0]T  in 
the  body  frame  while  the  target  docking  point  is  located  [1,  0,  0]  T  in  the  target  body 
frame.  The  initial  values  of  mentioned  and  the  remaining  states  for  computer  simulations 
discussed  in  the  following  sections  are  presented  in  Table  1.  The  scenario  also  assumes 
m  =  1  kg  ,  Q  =  0.005  rad/s ,  I '  =1  =  I3x3  (where  I3x3  is  the  identity  matrix),  r2  =1.5m  and 
the  maximum  bound  on  the  final  time  of  10  seconds.  The  choice  of  an  Q  value  was 
made  arbitrarily. 


State 

Initial 
condition 
(m  and  m/s) 

State 

Initial 

condition 

(rad/s) 

State 

Initial 

condition 

(quaternion) 

State 

Initial 

condition 

(quaternion) 

X\ 

0 

Xl 

0 

Xl3 

0 

X\1 

0 

X2 

5 

Xg 

0 

X\4 

0 

Xlg 

0 

*3 

0 

Xg 

0 

X\5 

0 

x\g 

0 

X4 

0 

Xio 

0 

X\6 

1 

X20 

1 

X5 

0 

xn 

0.25 

X6 

0 

Xl2 

0.25 

Table  1 .  The  initial  values  of  the  states 


2.  Solving  the  Problem  with  the  Gauss  Pseudospectral  Optimization 
Solver  (GPOPS) 

The  optimal  control  problems  posted  in  Section  A  were  first  solved  using  the 
Gauss  Pseudospectral  Optimization  Solver  (GPOPS)  (Rao  et  al.  2009).  The  GPOPS  is  an 
open  source  code  that  implements  the  Gauss  Pseudospectral  Method  for  solving  optimal 
control  problems.  As  many  as  150  internal  nodes  were  chosen  for  the  solution.  Usually, 
to  speed  up  the  numerical  procedure,  not  more  than  about  60  nodes  are  being  used 
(Yakimenko,  Xu  and  Basset  2008).  The  initial  conditions  were  chosen  based  on  Table  1 
and  the  final  conditions  were  based  on  matching  position  and  velocity  of  the  terminal 
point. 
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3.  Verification  of  the  GPOPS  Solution  with  Minimum  Principle 

Differential  equations  for  the  states  (Equations  (1),  (3),  (5),  (9))  and  costates 
(Equations  (22)-(25)),  developed  in  the  previous  sections,  were  programmed  into  the 
Mathworks’  Simulink  blocks  as  shown  in  Figure  4. 


Figure  4.  Simulink  model  for  integrating  the  states  and  costates  using  the  optimal 

switching  function  for  controls. 


The  Optimal  Control  Law  block  (on  the  left  in  Figure  4),  implements  the  control 
based  on  the  switching  conditions  Equation  (19)  or  Equation  (21)  developed  in  Chapter 
II.B.l.  The  upper  and  middle  blocks  integrate  the  states  and  costates,  respectively,  while 
block  on  the  bottom  calculates  the  Hamiltonian.  The  initial  conditions  for  the  20  states 
are  taken  from  Table  1,  while  the  20  initial  conditions  for  costates  along  with  the  final 
time  tf  are  guessed  on. 

The  far  right  block  calculates  the  terminal  conditions,  e  and  L, .  The  outputs  are 
sent  to  a  cost  function  ft  constructed  as  follows: 


/?  =  ||e(t/)||  +  ||^/)||  +  //(t/)2 


(47) 
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The  MATLAB  unconstrained  optimization  function  fminunc  (exploiting  Quasi-Newton 
method)  was  used  to  solve  for  the  TPBVP  in  the  following  fonn: 

X=fminunc ( @cost_function, xO, options, lambda, time) . 

This  function  attempts  to  find  a  minimum  of  a  scalar  function  of  several  variables, 
cost_function,  starting  at  an  initial  estimate  xO.  Here,  cost_function  calls  the 
Simulink  model  (Figure  4)  passes  the  vector  of  initial  guesses 

X  =  (48) 

and  computes  the  cost  function  p.  The  options  parameter  of  fminunc  is  varied 
based  on  the  fidelity  of  solution  required.  By  default,  the  termination  tolerance  the 
function  value  and  on  vector  of  varied  parameters  is  set  to  10"6. 

For  the  case  of  path  constraints  placed  on  the  state  variables  as  in  Equation  (31), 
20  additional  parameters  are  needed  to  define  the  costates  values  at  the  time,  td,  where 

they  become  discontinuous  (Bryson  and  Ho  1975).  For  this  case,  a  reset  is  included 
along  with  the  integrator  that  will  reset  the  costates  to  these  new  parameters  if  the  path 
constraint  is  contacted.  In  this  case,  the  X  value  is  augmented  to  include  initial  guesses  of 
the  costate  values  for  the  time  td  (Bryson  and  Ho  1975)  and  is 

X  =  [d1(t0),...,T20(t0),/i1(Q),...,i20(Q),t/]T  (49) 

The  vector  X  is  only  augmented  with  enough  costate  reset  values  as  suspected 
path  constraint  contact  points  deduced  from  the  GPOPS  solution.  For  example,  the 
GPOPS  solution  suggests  only  one  contact  with  the  path  constraint,  therefore,  we  only 
augment  the  X  vector  with  one  set  of  A (td)  values. 

To  implement  the  state  path  constraint  of  the  form  h(\)  >  0 ,  where  h  does  not 
depend  on  u,  a  penalty  function  P ,  was  associated  with  the  violation  of  the  constraint  that 
took  the  form  of: 
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(50) 


with 


[ 2xtx4  +  2x2x5  +  2x3x6  ,  h  <  0 
{  0,  h>  0 


P  =  J  pdt 


(51) 


and 


/3  =  \e{tf)\  +  \atf)\  +  H(tf)2  +P  .  (52) 

instead  of  Equation  (47).  This  increases  the  cost  function  if  the  vehicle  is  on  the 
constraint  boundary  and  not  meeting  the  tangency  conditions  stated  in  Chapter  II.B.2. 
Otherwise,  if  the  vehicle  is  not  on  the  constraint  or  meets  the  tangency  conditions  while 
on  the  constraint,  there  is  no  penalty  associated  with  the  cost.  A  penalty  function  of  this 
type  is  appealing  because  it  has  a  smooth  transition  from  solutions  where  constraints  are 
not  violated  to  solutions  where  they  are  violated.  Note,  that  if  the  vehicle  is  on  the 
boundary  and  meets  the  tangency  conditions,  it  will  not  cross  the  boundary. 

D.  OBTAINING  AND  ANALYZING  MINIMUM  QUADRATIC-CONTROL 
(MINIMUM  ENERGY)  SOLUTION 

The  results  of  the  GPOPS  solution  and  then  their  verification  with  the  Minimum 
Principle  solution  are  obtained  as  discussed  in  the  previous  section. 

1,  Minimum  Energy  Solution  with  GPOPS 

For  the  minimum  energy  rendezvous  scenario  set  in  Chapter  II. A,  the  GPOPS 
optimization  package  yielded  the  solution  shown  in  Figure  5.  This  solution,  minimizing 
the  quadratic-control  cost  returned  a  value  of  J  =  0.1133.  Figure  6  shows  the  planar 
views  of  the  solution.  The  overlaid  sphere  is  centered  a  the  target  RSO  and  has  a  radius 
equal  to  that  of  the  distance  the  docking  point  of  the  RSO  is  offset  from  its  center  of 
gravity.  The  final  maneuver  time  is  calculated  to  be  10  seconds,  the  upper  bound  on  the 
final  allowable  time  for  this  scenario  (without  this  limit  the  optimal  solution  would  yield 
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an  infinite  final  time).  It  should  be  noted  that  although  position  and  velocities  of  the 
docking  points  coincide,  the  orientation  of  the  spacecraft  do  not  match  since  this 
condition  was  not  set  as  a  constraint.  In  essence,  a  spacecraft  can  be  translating  at  a 
different  rate  than  the  docking  point  because  the  docking  point  is  offset  from  the 
spacecraft  CM  and  rotating  with  respect  to  the  orbital  frame. 


Figure  5.  Minimum  energy  solution  (GPOPS):  The  3D  view  of  the  optimal  trajectory. 

A  close-up  of  the  final  position  is  shown  in  the  exploded  view. 

The  progression  of  the  RSO  docking  point  position  carves  out  a  circle  in  the  plane 
perpendicular  to  the  angular  velocity  vector  of  the  RSO  because  the  RSO  is  assumed  to 
have  an  identity  inertia  matrix,  an  assumption  that  is  relaxed  in  the  next  chapter. 
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Figure  6.  Minimum  energy  solution  (GPOPS):  the  planar  views  of  the  optimal 

trajectory. 


Figure  7  shows  a  plot  of  the  resulting  controls  time  histories  (solid  lines)  as  well 
as  the  associated  costate  time  histories  (without  units)  that  were  used  to  synthesize  the 
optimal  control  based  on  Equation  (15)  (dashed  lines).  Figure  8  shows  the  time  history  of 
the  difference  in  position  and  velocity  of  the  chaser  and  RSO  docking  point  in  the  XYZ 
orbital  frame,  e(t),  as  defined  by  Equation  (12).  The  time  history  of  the  terminal 
conditions  due  to  transversality,  i '(t)  ,  are  shown  in  Figure  9.  Obviously,  all  the  values  in 
Figures  8  and  9  do  approach  zero  as  they  are  supposed  to  do.  Table  2  summarizes  the 
results  of  optimization  in  terms  of  the  values  of  varied  parameters,  initial  value  of  the 
costates  and  Table  3  lists  the  terminal  values  of  e(t),  £(t)  and  H(t).  Since  for 
numerical  solutions  the  final  value  of  the  Flamiltonian  does  not  tell  a  full  story,  its 
complete  time  history  is  presented  in  Figure  10. 
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Figure  7.  Minimum  energy  solution  (GPOPS):  control  and  associated  costate  history. 


Figure  8.  Minimum  energy  solution  (GPOPS):  time  history  of  discrepancies  in  the 

position  and  velocity  of  the  chase  and  RSO  docking  points. 
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Figure  9.  Minimum  energy  solution  (GPOPS):  history  of  the  transversality  conditions. 


Costate 

Initial  Condition 

Costate 

Initial  Condition 

Costate 

Initial  Condition 

X\ 

0.0202713560903818 

As 

0.0235328656506421 

Al5 

-0 . 00357471284206 

0.0415555942048028 

A9 

-0 . 058080907113755 

Al6 

1.060304309630e-10 

A3 

0.00694445082087804 

AlO 

-0.358237973748141 

Al7 

-0 .0681335012622741 

a4 

0.0816306658067988 

An 

0.556788661928892 

^18 

0 . 0053037488085563 

A5 

0 . 244702539834469 

A12 

0.410549229662363 

A19 

0.202974770150912 

^6 

0.0166838774111 

Al3 

1 ,41849853483e-05 

A- 20 

-3.422927080667e-10 

A? 

-0.004695462475917 

Al4 

-2 . 046214161772e-05 

10 

Table  2.  The  initial  values  of  costates  and  final  time  as  defined  by  GPOPS. 


endpoint 

Resulting  Value 

endpoint 

Resulting  Value 

endpoint 

Resulting  Value 

<?i  (tf) 

2.1  e-10 

cm 

-9.9  e-09 

CM 

7.7  e-08 

e2(G 

-9.4  e-12 

m 

-1.6  e-09 

Ciofi/) 

5.0  e-08 

?M 

-2.8  e-11 

cm 

1.2  e-11 

Cn(tf) 

1.4  e-07 

<?M 

5.1  e-11 

cm 

-4.5  e-09 

Cn(tf) 

6.5  e-08 

e5(f) 

1.7  e-10 

cm 

2.6  e-09 

CiM 

3.3  e-07 

1.4  e-10 

cm 

6.4  e-07 

C\M 

-5.8  e-07 

Cl  (0 

3.8  e-07 

CM 

-6.4  e-09 

H(tf) 

-0 . 009711 

Table  3.  Value  of  terminal  and  transversality  conditions  at  the  final  time. 
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Figure  10.  Minimum  energy  solution  (GPOPS):  history  of  the  Hamiltonian. 

It  should  be  noted  that  the  achieved  terminal  tolerance  of  the  order  of  1 0‘7- 10' 11 
does  not  necessarily  tell  about  the  quality  of  the  solution.  The  reason  for  this  is  that  the 
parameters  of  the  trajectory  are  only  being  computed  (optimality  conditions  enforced)  at 
150  nodes.  Another  important  issue  worth  mentioning  here  is  that  it  took  11,869.27 
seconds  (a  little  over  3  hours)  to  produce  this  150-node  solution  of  a  lOs-long  trajectory 
on  a  2.33GHz  Dell  Precision  M90  desktop  computer  with  an  Intel  T7600  processor  and 
1Gb  of  RAM.  The  initial  guess  for  the  solution  (required  as  an  input  to  GPOPS)  consisted 
of  two  terminal  points,  one  at  the  initial  time  and  one  at  the  final  time.  The  guess  for  the 
initial  states  corresponded  to  the  initial  conditions,  while  the  guess  for  the  final  states 
consisted  of  zeros  for  the  first  12  states,  and  the  value  of  [0,0,0, 1]  T  for  the  states 
corresponding  to  the  quaternions.  The  guess  for  the  control  history  was  0  at  initial  and 
final  time  for  all  controls.  Less  accurate  solution,  with  the  lower  number  of  nodes,  would 
obviously  require  less  computational  resources  (Yakimenko,  Xu  and  Basset  2008).  For 
instance,  a  25-node  solution,  for  this  case,  only  required  219.43  seconds  (less  than  4 
minutes)  on  the  same  computer.  The  point  is  that  the  solution  seems  to  be  feasible  and 
realizable  in  practice  (look  at  the  smooth  controls  in  Figure  7)  but  can  only  be  used  for 
off-line  computations,  i.e.,  in  the  open-loop  guidance  schemes.  The  jump  in  Hamiltonian 
value  in  Figure  10  and  other  Hamiltonian  histories  to  follow  occurs  when  path  constraint 
from  Equation  (32)  is  enforced  as  in  Equation  (36). 

The  next  section  addresses  validation  of  the  GPOPS  solution,  showing  that  it  is 
indeed  optimal. 
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2.  Verification  of  the  GPOPS  Solution  with  Minimum  Principle 

As  discussed  in  Chapter  II. C. 3,  the  initial  guesses  provided  by  the  GPOPS 
solution  (Table  2)  were  used  to  run  an  optimization  procedure  exploiting  the  optimal 
controls  synthesized  using  the  Minimum  Principle.  The  Quasi-Newton  method  based 
optimization  routine  employing  forward  shooting  and  integration  of  equations  of  motion 
using  Bogacki-Shampine  ordinary  differential  equation  3  (ODE3)  solver  with  a  fixed  step 
of  10'  .  This  approach  results  in  a  solution  that  has  as  many  as  10,000  points  (as  opposed 
to  just  150  nodes  as  in  GPOPS).  Of  course  it  comes  with  the  cost.  Even  with  the  perfect 
initial  guesses  for  all  varied  parameters  it  takes  many  hours  for  the  optimization  process 
to  converge  (compared  to  several  minutes  with  no  good  initial  guess  for  GPOPS).  The 
MP  solution,  returned  a  value  of  the  PI  ,7=0.1 1341 185  (compare  to  J  =  0.1133  for  the 
GPOPS  solution),  cost  function  p  =0.001768,  and  P  =  1.646  e-06. 

As  expected,  solving  the  same  problem,  as  in  Chapter  II.D.l,  (a  quadratic-control 
case  with  path  constraints)  produces  the  results,  the  optimal  trajectory,  controls  and  time 
histories  of  all  parameters,  which  are  pretty  close  to  those  produced  by  GPOPS.  The 
control  time  histories  and  those  states  and  costates  that  were  not  shown  for  the  GPOCS 
solution  are  shown.  Specifically,  Figure  1 1  shows  the  controls  time  histories, 
discontinuous  at  running  into  a  path  constraint  at  td  =9.6425 ,  Figures  12-14  show  the 
states  and  corresponding  costates  of  the  chaser  and  target  RSO.  The  values  of  the  initial 
costates  as  suggested  by  MP  are  shown  in  Table  4,  and  the  tenninal  values  e(tf) ,  C,(tf) 

and  77(C)  are  shown  in  Table  5.  Figure  15  presents  the  time  history  of  the  Hamiltonian. 
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Figure  1 1 .  Minimum  energy  solution  (MP):  Control  time  histories  and  associated  costates. 
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Figure  12.  Minimum  energy  solution  (MP):  state  and  costate  time-histories  for 

translational  variables  of  the  chaser  vehicle. 
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Figure  13. 
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Minimum  energy  solution  (MP):  state  and  costate  histories  for  the  defining 
angular  parameters  of  the  chaser  vehicle. 
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Minimum  energy  solution  (MP):  state  and  costate  time  histories  for  the  RSO. 
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Figure  14. 
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Table  4.  The  initial  values  of  costates  and  t/as  defined  by  MP. 


Costate 

Resulting  Value 

Costate 

Resulting  Value 

Costate 

Resulting  Value 

e\ (tf) 

0.0001 

m ) 

0.0001 

cm 

0.0001 

e2(tf ) 

-0.0001 

cm 

-0.0001 

Cio(h) 

-0.0001 

e3(h) 

-0.0006 

cm 

-0.0006 

Cn(tf) 

-0.0006 

e4(G 

-0.0033 

cm 

-0.0033 

Cn(tf) 

-0.0033 

es(tf) 

0.0012 

cm 

0.0012 

Cn(tf) 

0.0012 

eb(tf) 

0.0013 

cm 

0.0013 

Cu(t/) 

0.0013 

m 

-0.0001 

cm 

-0.0001 

H(tf) 

-0.0291 

Table  5.  The  value  of  terminal  and  transversality  conditions  at  the  final  time. 
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Figure  15.  Hamiltonian  for  the  forward  shooting  minimum-control  solution. 


Note,  some  of  the  numbers  in  Table  4  appear  in  the  bold  font.  These,  are  the  only 
numbers  that  are  different  from  the  GPOPS  solution  (Table  4).  Therefore,  it  took  hours 
and  hours  to  correct  initial  values  for  just  a  few  costates,  primarily  four  of  them,  Tn,  Tm, 
/li6,  and  a2o.  Obviously,  it  is  due  to  the  fact  that  because  if  integration  the  solution  is  quite 
sensitive  to  the  initial  values  of  varied  parameters.  Also,  the  accuracy  of  GPOPS  solution 
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(Table  3)  cannot  possibly  be  matched.  However,  the  point  of  the  entire  exercise  was  to 
prove  that  GPOPS  does  provide  a  solution  that  is  very  close  to  the  truly  optimal  one,  and 
it  was  proven. 

E.  OBTAINING  AND  ANALYZING  THE  MINIMUM  TIME  SOLUTION 
1.  Minimum  Time  Solution  with  GPOPS 

The  3D  trajectory  and  2D  projections  of  the  trajectories  as  a  result  of  the  GPOPS 
solution  are  shown  in  Figures  16  and  17,  respectively.  As  seen,  they  are  quite  different 
from  those  shown  in  Figures  5  and  6.  Also  notice  that,  although  the  positions  and 
velocities  of  the  angular  velocities  match,  the  vehicle  orientations  do  not.  The  resulting 
control  history  is  shown  in  Figure  18.  The  fz  control  (of  the  translational  motion  in  the  z 

direction)  turns  out  to  be  highly  oscillating.  The  endpoint  conditions,  taken  from  the  state 
constraints  on  the  endpoint  as  well  as  the  transversality  conditions  are  shown  in  Figure  19 
and  20.  The  values  of  the  initial  costates  as  suggested  by  GPOPS  are  shown  in  Table  6 
and  the  terminal  conditions  of  the  boundary  equations  are  shown  in  Table  7.  The  time 
history  for  the  Hamiltonian  is  presented  in  Figure  21. 
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x  (radial)  (m) 


Figure  16.  Minimum  time  solution  (GPOPS):  The  3D  view  of  the  optimal  trajectory.  A 

close-up  of  the  final  position  is  shown  in  the  exploded  view. 
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(jS/w)  (jS/w)  \  (jS/ui)  z; 


Figure  17.  Minimum  time  solution  (GPOPS):  the  2D  plots  of  optimal  rendezvous 

trajectory. 


Figure  18.  Minimum  time  solution  (GPOPS):  the  control  and  associated  costate  history. 
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Figure  19.  Minimum  time  solution  (GPOPS):  time  history  of  discrepancies  in  the 
position  and  velocity  of  the  chaser  and  RSO  docking  point. 


Time  (s)  Time  (s) 


Figure  20.  Minimum  time  solution  (GPOPS):  history  of  transversality  conditions. 
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2, 

Initial  Condition 

2, 

Initial  Condition 

2, 

Initial  Condition 

0.000834092040702561 

>- 

OO 

-0 . 00781391345181247 

2-15 

-0.0670307214682799 

^2 

0 .418560914087817 

Ag 

-0 . 00529954719645793 

2-16 

-2.47014407229089 

^3 

7 . 46232139969827e-08 

2io 

-0 .0966015054618492 

2-17 

-0.365359323439127 

A4 

-0 . 00634589252204042 

In 

0 . 017598938378927 

OO 

0.0400467496481145 

^5 

0 . 901334242652277 

2 12 

-0.719562963054045 

2-19 

-0.562241844247041 

/t.6 

1 . 04649495108653e-07 

2l3 

-0 . 0420682846405162 

220 

-1 . 6978671241818 

A? 

-0 . 0133948065189441 

2l4 

-0 . 0394802339637517 

tf 

3.4237 

Table  6.  The  initial  values  of  costates  and  as  defined  by  GPOPS. 


endpoints 

Resulting  Value 

endpoints 

Resulting  Value 

endpoints 

Resulting  Value 

ei(tf) 

-3.1  e-08 

cm 

3.3  e-07 

cm 

-1.9  e-06 

ei(tf) 

-3.8  e-08 

cm 

-6.9  e-08 

Cioth) 

3.8  e-06 

e3(f) 

4 . 9  e-08 

cm 

-5.2  e-14 

cm 

3.2  e-07 

e4{tf) 

-1.8  e-07 

cm 

1.0  e-07 

Cn(tf) 

-5.1  e-07 

<?5  (tf) 

5.1  e-07 

cm 

1.0  e-07 

cm 

-2.0  e-07 

e6{tf) 

-1.7  e-08 

cm 

-3.5  e-07 

cm 

1.1  e-07 

cm 

-3.3  e-10 

cm 

9  7  0— 0  6 

H(tf ) 

-0.2801 

Table  7.  Value  of  terminal  conditions  at  the  final  time  as  calculated  by  GPOPS. 


Figure  2 1 .  Minimum  time  solution  (GPOPS):  history  of  the  Hamiltonian. 


Closer  inspection  of  the  GPOPS  solution  for  in  Figure  18  (it  jumps  back  and 
forth?  around  zero  value)  suggests  a  presence  of  a  singular  control.  This  is  indirectly 
confirmed  by  oscillations  of  the  Hamiltonian  as  well  (Yakimenko,  Xu  and  Basset  2008). 
Hence,  the  fz  control  is  not  optimal  and  infeasible. 

With  150  nodes,  the  computational  time  to  arrive  at  the  solution  shown  above 
(3.4237  second  maneuver)  was  8,929.55  seconds  (~2.5  hours).  The  initial  guess  for  the 
solution  was  the  same  as  the  one  used  for  the  minimum  energy  case.  For  comparison, 

with  25  nodes  the  required  computational  time  can  be  brought  down  to  100.77  seconds 
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(<2  min).  However,  as  seen  from  Figure  22,  the  controls  in  this  case  are  even  less 
optimal  than  in  Figure  18,  and  of  course  the  Hamiltonian  (Figure  23)  looks  worse  than 
that  of  Figure  21. 


Figure  22.  Minimum  time  solution  (GPOPS):  The  associated  costate  history  resulting 

from  a  25  node  solution. 


Figure  23.  Minimum  time  solution  (GPOPS):  Hamiltonian  for  the  associated  25-node 

solution. 

2.  Verification  of  the  GPOPS  Solution  with  Minimum  Principle 

The  minimum  time  rendezvous  problem  is  approached  in  a  similar  fashion  of  the 
minimum-control  one.  First,  the  problem  is  investigated  using  a  ft  based  on  Equation 
(52).  A  major  difference  that  arises  compared  to  the  minimum-control  solution  is  the 
existence  of  a  singular  control  in  M3,  which  controls  acceleration  in  the  z  orbital  direction. 
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As  discussed  in  Chapter  II.B.3,  a  singular  control  exists  and  needs  special  consideration 
for  numerical  implementation.  The  optimal  control  structure  was  assumed  to  have  the 
following  form: 


w3 


-1,  0<t<tx 
0,  tx<t  <t2 
1,  t2<t  <t 3 
<  0,  t3<t  <t4 
-1,  t4<t  <t5 
0,  t5<t<  t6 

1,  t6<t<tf 


(53) 


Accordingly,  the  vector  of  varied  parameters  X  was  augmented  with  the  switching 
times  ; 


X  —  \\(tf ),•••, A2Q(tf),Ax(td),...,A2Q(td),tx,...,t6,tj-]  .  (54) 

Having  instances  tx,...,t6  as  varied  parameters  implied  that  the  search  was  made 
among  multiple  control  profdes  including  traditional  bang-bang  control,  like  w3min  -u3max 
(when  tx=t2  and  /,.  =  tf ,  i  =  3,..., 6),  or  w3max-w3min  (when  tt  =t2=  0 ,  t3=t4  and 
t5  =t6=  lf  ,  i  -  3,..., 6)  allows  calculation  of  a  control  history  even  in  the  presence  of  a 
singular  arc. 

The  resulting  trajectory  is  shown  in  Figures  24  and  25.  The  corresponding  optimal 
controls  profdes  are  presented  in  Figure  26.  As  seen,  it  was  w3min  -  0-w3max  -w3min  profde 
that  was  found  to  be  optimal  for  the  w3  control  (i.e.,  t5=t6  =  tf  ).  The  time  histories  for 
the  remaining  controls,  ux,u2,u4,u5,u6,  match  those  for  the  GPOPS  solution  (Figure  15). 
Figures  27-29  show  the  state  costates  time  histories  for  the  chaser  and  target  RSO, 
respectively.  The  values  of  the  initial  costates  as  suggested  by  MP  are  shown  in  Table  8 
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x  (radial)  (m) 


and  the  error  in  satisfying  boundary  conditions  are  shown  in  Table  9  (the  cost  /?  =  0.002, 
P=0.00124  Figure  30  presents  the  Flamiltonian  (compare  it  with  that  of  Figure  21  and 
22). 


1  i'i 

0.5 -■ 

o  — 

0.5  -J  i 


'Target  docking  point 
■  Chaser  CG 
•  Chaser  docking  point 


y  (in-track)  (m) 


-u.o 

z  (cross-track)  (m) 


Figure  24.  Minimum  time  solution  (MP):  3D  optimal  rendezvous  trajectory. 
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(jS/uj)  *1  (jS/m)  \  (jS/ui)  ZJ 


0  1  2  3  4  5  -1  -0.5  0  0.5 

y(m)  z  (m) 


Figure  25.  Minimum  time  solution  (MP):  The  2D  projections  of  the  optimal  rendezvous 

trajectory. 


Time  (s)  Time  (s) 


Figure  26.  Minimum  time  solution  (MP):  control  histories  and  associated  switching 

conditions. 
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Figure  27. 


Figure  28. 
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Minimum  time  solution  (MP):  State  and  costate  histories  for  the  translational 
variables  of  the  chaser  vehicle. 


Time  (s)  Time  (s) 

Minimum  time  solution  (MP):  state  and  costate  histories  for  the  defining 
angular  parameters  of  the  chaser  vehicle. 
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Figure  29.  Minimum  time  solution  (MP):  state  and  costate  time  histories  for  the  RSO. 


2, 

Initial  Condition 

2, 

Initial  Condition 

2, 

Initial  Condition 

k 

0.000813270815835S08 

^8 

-0.0079043944784984 

An 

-0.067030975848676 

fa 

0.41856021740054 

-0.0052514744594656 

Al6 

-2.47013927846309 

k 

0 

•^10 

-0.0966015055875181 

An 

-0.365359450504476 

A4 

-0.00632550328119496 

^11 

0.017598938703864 

An 

0.0400468901069352 

A5 

0.90134410424794 

A\2 

-0.719562962786195 

A\9 

-0.562241768137901 

Ae 

0 

An 

-0.042050296747438 

^20 

-1 . 6978671241818 

I7 

-0.0135549835134852 

A\4 

-0.0395088309550527 

h 

3 . 432 

Table  8.  The  initial  values  of  costates  and  tf  as  defined  by  MP. 


e, 

Resulting  Value 

e, 

Resulting  Value 

e, 

Resulting  Value 

ei  Uf) 

0.016 

cm 

0.0028 

cm 

0.0096 

-0 . 014 

cm 

-0.0018 

Cioth) 

0.0061 

e3fi/) 

0.0016 

m > 

0.00013 

Cnth) 

0.00048 

£4  (tf) 

-0 . 017 

cm 

-0.00059 

Cu(tf) 

-0.00079 

es(tf ) 

0.029 

cm 

-0 . 00077 

Cn(tf) 

0.0010 

ee(tf) 

0.0080 

cm 

0.0066 

Cu(tf) 

6 . 6  e-05 

cm 

1 . 9  e-05 

cm 

-0.0063 

H(tf) 

0 . 15 

Table  9.  Value  of  terminal  and  transversality  conditions  at  the  final  time. 
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Figure  30.  Minimum  time  solution  (MP):  Flamiltonian  for  the  forward  shooting  minimum 

time  solution. 

Again,  the  bold  numbers  in  Table  8  indicate  the  differences  as  compared  to 
GPOPS  quasi-optimal  solution.  As  opposed  to  the  minimum-control  case,  when  truly 
optimal  solution  happened  to  have  almost  the  same  values  of  the  initial  costates,  in  the 
minimum  time  case  implying  a  singular  control,  the  optimal  solution  involved  more 
variations  from  the  solution  provided  by  GPOPS.  Of  course,  as  in  the  minimum-control 
case  presented  in  Chapter  II. D. 2,  the  indirect-method-based  optimization  was  run 
recursively  multiple  times  for  several  days  in  order  to  converge.  Nevertheless,  the  errors 
in  satisfying  the  terminal  conditions  (Table  9)  were  of  several  orders  higher  than  that  of 
claimed  by  the  GPOPS  solution  (Table  7).  The  next  section  provides  some  more  details 
on  this  issue. 

3.  Propagation  of  the  GPOPS  Solution 

Once  again,  for  the  solutions  provided  by  GPOPS,  vehicles  dynamics  are  only 
satisfied  in  a  limited  number  of  nodes.  That  is  why  the  error  in  meeting  all  constraints  is 
so  negligibly  small  as  compared  to  solutions  provided  by  the  methods  involving 
integration  of  equations  of  motion  (shooting).  Flowever,  if  we  integrate  the  GPOPS 
controls,  we  will  end  up  with  about  the  same  (lower)  accuracy  as  the  shooting  approach. 

To  illustrate  it,  let  us,  for  example,  take  the  controls  for  the  minimum  time 
solution  (including  a  “weird”  control  for  /))  and  integrate  it.  The  result  of  integrating  the 
equations  of  motion  derived  in  Chapters  II.  A  and  control  switches  in  Chapter  II.B  with  a 
fixed  time  step  of  0.0001  seconds  is  shown  in  Figures  31  and  32  (the  zero-order  hold  of 
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the  last  control  inputs  was  used).  As  seen,  the  trajectories  are  practically  the  same  as  in 
Figures  16  and  17,  but  the  endpoint  discrepancies,  summarized  in  Table  10,  obviously 
grew  up.  The  endpoint  conditions  of  the  state  variables  are  shown  in  Table  10.  Note,  that 
the  endpoint  conditions  of  the  transversality  conditions  are  not  known  because  the 
costates  are  not  propagated. 


Figure  31. 


Propagated  trajectory  using  the  minimum  time  control  history  from  GPOPS. 
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y(m)  z  (m) 


Figure  32.  Propagated  trajectory  using  the  minimum  time  control  history  from  GPOPS. 


e 

Resulting  Value 

e 

Resulting  Value 

e 

Resulting  Value 

-0.00027 

Clitf) 

N/A 

cm 

N/A 

^2  (fr) 

-9 . 3e-05 

cm 

N/A 

Cio(fr) 

N/A 

^3  (tf) 

0.00037 

cm 

N/A 

cm 

N/A 

e4(fr) 

-0.00092 

cm 

N/A 

cm 

N/A 

<?5  (tf) 

-0.00035 

cm 

N/A 

cm 

N/A 

0.00026 

cm 

N/A 

Cu{tf) 

N/A 

cm 

N/A 

cm 

N/A 

mtf) 

N/A 

Table  10.  Value  of  terminal  conditions  at  the  final  time. 
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III.  OPTIMIZATION  OF  A  SPACECRAFT  MANEUVER  FOR 
CLOSE  APPROACH  AND  DOCKING  WITH  A  TUMBLING 

OBJECT 


This  chapter  takes  the  study  further  by  considering  translational  actuators  fixed  in 
the  chaser  spacecraft  body  frame  (better  representing  an  actual  system),  adding  the 
requirement  to  match  attitude  and  angular  velocity  of  the  chaser  and  RSO  in  preparation 
for  a  docking  procedure  (along  with  the  requirement  of  matching  position  and  velocity  of 
chaser  and  RSO  docking  points),  and  considering  one  more  PI  (minimum  fuel).  The 
optimal  control  problem  posed  in  Chapter  II  is,  therefore,  reformulated  to  take  into 
consideration  the  effects  on  dynamics  and  the  extra  constraints. 


A.  RENDEZVOUS  MODELING  AND  OPTIMIZATION  PROBLEM 
FORMULATION 


Using  these  controls,  we  would  like  to  bring  the  two  spacecraft  from  some  initial 
conditions,  given  by  20  initial  values  of  states  xi(t0) ,  /  =  l,...,20,  to  a  docking-enabling 
condition  described  by  matching  the  chaser’s  and  RSO’s  docking  station  final  positions 
and  velocity  vectors  as  well  as  matching  the  orientation  and  angular  velocity  of  the  two 

vehicles.  Defining  °T  R  as  the  rotation  matrix  to  convert  from  the  body  frame  of  the  target 
to  the  orbital  frame,  we  can  express  the  desired  terminal  conditions.  These  conditions  of 
matching  spacecraft  docking  positions  and  velocities,  in  the  matrix  form,  are  identical  to 
the  conditions  shown  in  Equation  (12).  Matching  orientation  and  angular  rates  leads  to 
seven  more  boundary  equations  at  tf. 


1 

1 

1 

_ 1 
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i _ 

T 

q2 

q2 

^8 

T 

q^ 
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I 

_ 1 
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and 


55 


(56) 


r 

r  c 

-1 

cox 

\  COx  ! 

ell 

„ T 

r 

co2 

— 

o2 

— 

e\2 

r  1 

c 

0)3  1 

co3  1 

_^13_ 

If  Equations  (55)  and  (56)  are  satisfied,  that  would  provide:  =  ° R  and  °(0<:  =°  to/ . 

Therefore,  the  conditions  in  Equation  (12)  can  be  rewritten  as: 


and 


~d\  -df 

X 

dl-dl 

- 

y 

= 

e2 

1 

^3 

1 

E— i  ^ 

^3 

_ i 

z 

_e3_ 
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-df 

X 

^4 

o  ,jt  _  On 

0)  x  rK 

dl 

-dc2 

- 

y 

= 

dl 

1 

^3 

1 

z 

_e6_ 

(57) 


(58) 


The  resulting  complete  set  of  desired  final  conditions  are  [e,  (tf), e13(//  )]T  =  0  . 

For  optimal  control  problem  in  this  chapter,  three  different  perfonnance  indices 
each  of  the  (Mayer)  form 

J=\fodt  (59) 

h 

with  the  running  cost 


/o=l 


(60) 


for  minimum  time, 
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(61) 


fo  -  ~^{U\  +U22  +  Uf  +  Uf  +U52  +U6 ~) 

for  minimum  quadratic-control  expenditure  (i.e.,  minimum  energy),  and 

fo  =  (|wi  |  +  |m2  |  +  |w3  |  +  |w4  |  +  |w5 1  +  |w6 1)  (62) 


for  minimum-control  are  considered.  Each  case  has  a  maximum  finite  value  chosen  for 
the  final  time  t f  but  is  still  considered  a  free  (varied)  parameter.  The  next  step  is  to 

convert  this  optimal  control  problem  to  a  TPBVP  using  MP  (Pontryagin  et  al.  1964; 
Bryson  and  Ho  1975) 

B.  MINIMUM  PRINCIPLE  FORMULATION 

This  section  deals  with  the  MP  formulation  and  synthesis  of  the  optimal  control. 
Using  the  same  state  vector  of  the  complete  system  as  shown  in  Equation  (10),  (the 
differential  equations  for  parameters  of  chaser  and  target  quaternions  are  standard  and 
can  be  found  in  Chapter  II),  and  the  control  vector  shown  in  Equation  (11),  we  write 
down  the  Hamiltonian: 

H  (k,  x,  u)  :=  f0  +  (k,  x)  +  juh  (63) 

where  operator  (...)  on  the  right-hand  side  denotes  a  scalar  product  of  two  vectors,  and 
k  e  tR  v’  is  a  costate  vector  where  its  differential  equations  are  to  be  defined  later  in  this 
section.  The  value  of  n  is  a  constant  and  is  dictated  by  the  same  constraint  as  in  Equation 
(36)  when  considering  path  constraints. 

Again,  as  in  Chapter  II.B.2,  this  constraint  is  that  the  CM  of  the  chaser  spacecraft 
must  remain  at  a  distance  larger  than  some  minimum  distance  (a  “keep  out”  sphere  with  a 
radius  r)  from  the  CM  of  the  RSO,  which  is  coincident  with  the  origin  of  the  orbit  frame. 
This  assures  that  the  chaser  vehicle  will  not  pass  through  the  target  vehicle  in  order  to 
reach  the  docking  position  shown  mathematically  in  Equation  (32). 

The  entire  Hamiltonian  is  expressed,  with  respect  to  the  state  vector  x  defined  by 
Equation  (10),  in  Equation  (64). 
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H  (k,  x,  u)  =  1  +  2,x4  +  A2xs  +  A3x6 

24(2fix5  +  3fi2x,  +  (x,26  +  x,23  -x,24  -x,25 W  +  2(x13x14  ~x{ixX6)n2  +  2(x13x15  +  x14xI6)w3) 


H - 

m 


+/l5(  2fix4  +  2 (x13x14  +  ) Wj  +  (x16  xu+xl4  X\s}u2  +  2(x14x15  x13x16)i/3) 

+26(-fi2x3  +2(x13x15  -x14x16)z/,  +  2(x14x15  +  x13x16)z*2  +  (x,26  -x2  -x,24  +  x2  )w3) 


v 

Ic  -Ic 

^  O')  1 


Ic-Ic  Ic-Ic 

+  ~  c  ~  (*8*9  +  «4  )^7  +  33Tr  ~  (*7*9  +  «5  H  +  ‘'r  "  (*8*7  +  M6  H 


lll 


*22 


+  II2  J33  (*n*i2Ho  +  /ii.r/"  (*io*i2Hi  + 
Al  ^22 


+2,3  — 

13  2 


IT-IT 

^7^(*n*10H: 


+2, 


14 


*9  -((*16)2  -(*13)2 -(*14)2  +(*l: 

-(*8  -2  (*14*15 -*13*16  )Q)  *15  +(*7  -2  (*13*15  +*14*16  )Q)*1, 
(*9  ((*16  )”  -(*13)2  -(*14)2  +(*11 

+  (*8  -2  (*i4*i5  -*13*16  )Q)  *i6  +(*v  -2  (*13*15  +*i4*i6  )Q)*i; 


fi  x 


fi  x 


+2,5  — 

15  2 


(*7  2  ( 


*13*15  +  *14*16 


)«)- 


+ 


(*8  2(x,4X,5  ,X13X16)Q^X13 +|x9  ^(x,6)  (*i3)  (*14)  +(*11 


fi  x, 


+26  — 
16  2 


+/l17- 


+4  2 


-(*7 -2( 


*13*15  +  *14*16 


)n) 


X- 


-(*8  2(x14x15  x13x16)Q)x14  |x9  ^(x16)  (*i3)  (*i4)  +  (* 

*12  ^(*20 )  —(*n)  —  (*is)  +(*] 


fix, 


fix, 


(*11  2(x18x19  xi7x2o  )  fi)  *i9  +  (*10  2(x17x19  +  .x18x20)fi)x2 


|*i2  ((*20)  (*17)  (*is)  +  (*11 


fi  X, 
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(xu  2(x18x19  x17x20)fi^x20  +(x10  2(x17x19  +x18x20)fi)x19 


+^  2 
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*17*19  +  *18*20 
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+ 


(*11  2  (*18*19  *17*2o)^)  *17  +|*12  ((*2o)  (*n)  (*18)  +(*11 


fi  X. 
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X,- 
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(64) 
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For  the  minimum  time  cost  problem  defined  by  Equation  (60),  the  part  of  the 
Hamiltonian  that  depends  on  the  controls,  the  switching  function,  turns  out  to  be: 

Ht('k,x,\x)  =  —(XJx+XJ  +A6f:)  +  ^A7u4  +-^\u5  +-^rA9u6  (65) 

m  hi  1  22  1 33 


where 


f,  =  ( cu  +  q? 2 - qh 2 ~(i32)u i  +  2 ( dUC2  ~ % q\ ) ■ u2  +  2 ( qfq'  +  ^ qh ) u3 
fy  =  2  (qciqC2  +  qU* ) w.  +  (q?  -  q?2  +  qf  ~  qf )  u2  +  Afif  -  q?q* )  u3  ■  (66) 

fz  =2(^1C^3C  ~q<2.q(A  )Ul  +2(^2^3C  +<7iC<74  )M2  +(^4^  “  ^2  "  + 


This  results  in  Equation  (65)  written  as: 


*  1  ,  .111 
H  (k,x,u)  =  — (switchjWj  +  switch,w,  +  switch,w, J  +  —  A1u4  +  —  \u5  +  —  A9u6  (67) 


m 


22 


33 


where  the  switching  functions  that  defined  the  control  above  are  given  by  the  relations: 

switch,  =  T4  (qf  +  qf  -  qf  -  qf  )  +  As2  (qcxq^  +  q^q^  )  +  \2  (qfq%  -  qc2qc4  ) 
switch,  =A42( qfq£  -  q*qc4  )  +  T5( qf  -  qf  +  qf  -  qf )  +  A6  2  ( q^q^  +  q^qc4  ) .  (68) 

switch,  =  A4  2  (q*q*  +  q7q4  )  +  As  2  (  q7q4  -  qfq4  )  +  A6  (qf  -  qf  -  qf  +  qf  ) 


Assuming  there  are  no  singular  arcs,  since  all  six  controls  enter  the  switching 
function  (Hamiltonian)  linearly,  the  optimal  control  that  minimizes  the  Hamiltonian  in 
the  case  of  minimum  time  is  the  bang-bang  control  defined  by: 


u\ 

U4 


!  1,  switch,  <  0 
[-1,  switch,  >0 

f  1,  A7<0 
|-1,  f>0 


1,  switch,  <  0 

<  iu 
[-1,  switch,  >0 

f  i,  K  <  0 

<  u, 

1-1,  ^>0 


;  1,  switch,  <  0 
[-1,  switch,  >0 

1  1,  f<0 

l-i,  ^>0 


(69) 


Likewise,  developing  the  Hamiltonian  for  the  minimum  quadratic-control  PI  based  on  the 
running  cost  (61),  the  switching  function  becomes: 
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(70) 


H  (X,,  x,  u)  —  —  (mj  +  u0  +  +  w4  +  W5  +  Wg  ^ 

+  — (4./l  +/l5/2  +^6/3)  +  -^^W4  +7^AM5  +7T^9M6 


22 


33 


Taking  into  account  Equation  (66),  the  resulting  optimal  control  that  minimizes  the 
Hamiltonian  for  the  minimum  quadratic-control  (minimum  energy)  PI  becomes: 


U;  =  i 


switch - 


m 


<-l 


switch 

m 

-1, 


and  u,  = 


'+3 


-1, 


/  =  1,2,3 


m 

switch; 
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>1 


A 


'i+3 


/ 


<-l 


-1<^<1  i  =  4,5,6,  7  =  /-3 


(71) 


A 


+3 


IL 


>1 


For  the  minimum  fuel  case,  corresponding  to  the  running  cost  in  Equation  (62),  the 
switching  function  can  be  written  as: 


H*(k,x,u)  =  (|  iq  |  + 1  u2 |  +  |w3 1  + 1  u4  I  +  \u5 1  +  |w6  1^ 


+  (^4/1  “*“'^'5 2/2  “*"^'6/3)"*"  rc  /^7W4  rc  AM5  rc  ^9W6 

m 


(72) 
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and,  therefore,  the  optimal  control  structure  becomes: 
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u,  =  < 


and  u,  = 


1, 

0, 

-1, 

1, 

0, 

-1, 


switch;  <  -m 

-m  <  switch;  <  m  i  =  1,2,3 
switch;  >  m 

4+3  <  '  4/ 

-/J<4+3</J  i  =  4,5,6,  7  =  1-3 

4+3  >45 


(73) 


The  differential  equations  for  costates  are  obtained  via  differentiating  the  Hamiltonian: 


SH 

Sx 


=  -)„ 


(74) 


which  yields: 


4  =-3A4D.  m  \  4  =0,4  =A6Ci2m  1 

=  ~4  +  24 Om  1 ,  4  =  -4  -  24 £lm~l  .  (75) 

4  =  ~4 


The  next  six  adjoint  equations,  corresponding  to  costates  7-12,  take  the  form  of 
Equations  (76)— (8 1): 


Ic -Ic  Ic  -Ic 

^  _  -'ll  -*33  3  v  1  1 22  -'ll 
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IC  ~IC  IC  -IC 

4  =  33/C  22  4x9  +  22/C  11  4X7  +  43  2  ^15  "  44  2  ^16  "  4s  2  ^13  +  4e  2  Qxi4 


Ic -Ic  Ic -Ic 

2  _  4 33  '22  1  ,  -'ll  4 33 

/L9  —  rC  ^^g 


c  4x7  43  ~  ^-14  +  44  ~  ^X|3  4s  ~  ^-16  +  4fi  ~  ^-xi5 


‘22 


I  -I  I  -I  1  1  1  1 

4o  =  Jr  —  4lX12  "l  Jr  42X11  _  4?  'J^X20  _  4s  'J^X19  +  49  J'^X18  +  4o  'J^X17 


‘22 


‘33 


/  -/  I  -I  1111 

4l  =  Jr  4oX12  ^  Jr  42X10  "*"47  J'^X19  _  4s  J'^X20  _  49  “f^i7  +  4o  J'^X18 


‘33 


(76) 

(77) 

(78) 

(79) 

(80) 
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f  -f  IT -IT  1  1  1  1 

A2  ~  Ao*n  "l  Tt  ~ Ai*io  —  A7  As  n  —  A9  ^'^*20  "*"  Ao  ~£^i9  ($1) 


ln 


22 


The  next  four  equations,  corresponding  to  costates  13-16,  are  of  the  form  of  Equations 
(82)— (85): 
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m 


f  A4(2xnul  +  2 xuu2  +  2x15w3)  +  /t5(2x14w1  -2 x13w2  -2x16w3) 
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1 
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■As  ^(_(x7  -2 (*13*15  +  *14*16)°)  +  2*16^*14  ~ 2x15Qx13  +  2xi4Qx,6  ) 
■A6^-(2xi6f2xi3  -(xg  -2(x14x15  -x13x16)Q)  +  2x15Qx14  -2x14Qx15) 


(83) 
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2,5 =— 
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^+A6(2x13u1  +  2x14w2  +  2x15m3) 
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m 
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-41(-(2x1Jnx14-2,1Jnx14+2x1!n,16) 


(85) 


and,  the  final  four  adjoint  equations  for  costates  17-20,  take  the  form  of  Equations  (86)- 
(89): 


2,7  =  -An  ^-(2x17Qx18  -2x20Qx19  -2x19Qx20) 


"  2,8  — ('T12  ~~  (*20  ~Xn  "^18  +-^19)^2  j 
"  2,9  ~(2x19Qx18  +(xn  —  2(xlgx19  — x17x2 
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As  - 
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-  4  ^-(2x17Qx18  -  2x18Qxp  +  2x19Qx20  ) 


-  Ao  ^-(+2*i720  +  2x182Q  -  (x12  -  ( 


2  _  2  _  2  2 

I  x20  x17  x18  -I-  x19 


)q)  +  2x192Q 


Ao  —  A7  2  (  2x20Qx18  2x17Qx19  +  (x10  2  (x17x19  +  x18x20 )  Q)  2x2or2x18j 

-  As  ^-(2*20^*17  +  2x20Qx17  +  (xn  -  2  (x18x19  -  XpX20 )  Cl)  -  2x18Qr19 ) 

-  Ao  ^(2*i72^  +  2x182Q  -  (x12  -  (x20  -  Xp  -  xf8  +  xf9 )  Q)  +  2x202n( 
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-  Ao  -(2xpQxi8  -2xpOx18  +2x19Ox20) 


(87) 


(88) 


(89) 


The  form  shown  in  Equations  (75)-(89)  is  different  than  presented  in  the  previous 
chapter,  again  because  of  the  body-mounted  thrusters. 

The  next  step  is  to  address  the  transversality  conditions,  which  with  account  to  the 
following  coupled  tenninal  variations,  followed  directly  from  Equations  (55)— (58).  The 
ten  scalar  equations,  eve2,...,e10,  represent  the  matching  of  translational  and  angular 

position  and  velocity  of  the  docking  points  in  the  three  components  of  the  orbital 
coordinate  frame  at  the  time  t  ,  based  on  Equations  (57)  and  (58).  Variable  E  represents  an 

endpoint  cost  that  is  set  to  zero  for  both  the  minimum  time  and  minimum-control 
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scenarios  (Yan  et  al.  2002)  and  are  equivalent  to  Equations  (26)  and  (27)  of  Chapter  II. 
The  first  six  equations  result  in  expressions  that  only  contain  the  parameters  t;  , 
i=l,. ..,6. 


f)  C  5  f)  t>2,  -^(C)  ^35  -^(fo)  C45  A'sif  f)  ^5’  ^(c)  ^6  (90) 


Realizing  the  specific  relationships,  such  as: 


SE 


fix7 


=  Kjf  +  V 7  and 


fi£ 


fix 


10 


(91) 


and  combining  terms,  seven  more  conditions  are  developed  that  must  be  satisfied,  along 
with  six  docking  states  for  position  and  velocity  given  by  Equations  (57)  and  (58)  and 
the  seven  angular  states  based  on  Equations  (55)  and  (56),  atr;. .  These  conditions  are 

shown  in  Equation  (92): 

(//  )  =  ^“7 tf  +  \otf  +  fl  (Xi/  ’  ^tf  )  +  flO  (X//  »  ^7/  )  =  ei4  (^/  )  =  9 
^56  (^/)  =  ^77/  +^(;+3)r/  +  +  f(i+3)(Xtf’^‘tf)  =  e(i+7)(^/)  =  fo  f°r  Z 

^7  (f/) = ^i?/- + +y/(xr/5  ) +y<;+4)(x?/7  ^i/) = e(i+4>(^/) = 9?  for  1 


=  7,8,9  and 
=  13,.  ..,16 


(92) 


that  only  depends  on  x;/  and  /.;/  .  This 

results  in  the  remaining  final  conditions,  C, ,  along  with  Equations  (55)— (5  8)  that  need  to 
be  satisfied  for  an  optimal  trajectory,  £  =  [£j (r, g7 (tf  )]T  =  0  . 

Other  parameters  for  the  scenario  are  m  =  1kg,  Q  =  0.005  rad/s ,  \T  =diag  ([3,1,2]), 
I1  =  l3x3  (where  l3x3  is  an  identity  matrix)  with  inertia  units  of  kg-m2  and  tf  set  to  a 

maximum  of  10  seconds.  This  scenario  initializes  the  chaser  CM  starting  at  a  distance  of 
5  m  from  the  target  CM  with  the  target  having  an  initial  angular  velocity  of  0.25rad/s  in 
both  y  and  z  body  axis.  The  body  coordinate  frames  of  each  spacecraft  and  the  orbit 
frame  are  assumed  to  be  coincident  with  the  inertial  frame  at  the  beginning  of  the 
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where  /(x,,,k(/.)  is  the  portion  of 


f  8EV 

dx 


simulation.  The  docking  point  of  the  chaser  spacecraft  is  located  at  0.25  m  in  the  negative 
x  body  frame  and  the  docking  position  of  the  target  is  located  at  lm  in  the  positive  x  body 
frame.  Therefore,  the  orientations  of  the  two  spacecraft  will  be  coincident  at  docking.  The 
resulting  values  of  the  initial  states  for  computer  simulations  discussed  in  the  following 
sections  are  taken  from  Table  1  in  the  previous  chapter. 

C.  SOLVING  THE  PROBLEM  NUMERICALLY  USING  A 

PSEUDOSPECTRAL  METHOD 

The  optimal  control  problem  posted  in  Chapter  III. A  was  again  solved  using 
GPOPS  (Huntington  and  Rao  2008).  This  software  solves  for  the  optimal  control,  state 
and  costate  history  based  on  a  given  PI  and  constraints.  As  opposed  to  the  150  node 
solutions  of  the  last  chapter,  200  nodes  were  chosen  for  the  solution  to  be  consistent  with 
previous  research  on  this  topic  (Ma  et  al.  2007). 

1.  Minimum  Time 

The  3D  trajectory  and  2D  projections  of  the  trajectories  are  shown  in  Figures  33 
and  34,  respectively.  In  these  figures,  the  solid  line  with  the  square  markers  shows  the 
trajectory  of  the  docking  point  aboard  the  target  RSO,  the  dotted-dashed  line  with  a 
circular  marker  represents  the  chaser  spacecraft’s  CM,  and  the  line  with  upside-down 
triangles  depicts  the  trajectory  of  the  docking  point  aboard  the  chaser  spacecraft.  The 
overlaid  sphere  is  centered  at  the  target  RSO’s  CM  and  has  a  radius  equal  to  that  of  the 
distance  the  docking  point  is  offset  from  its  CM. 

The  resulting  control  history  is  shown  in  Figure  35.  Figures  35-38  show  the  state 
and  costate  histories.  Figures  39  and  40  show  the  final  difference  in  endpoint  conditions 
approaching  zero  and  Figure  41  shows  the  time  history  of  the  transversality  conditions, 
Equation  (92),  as  they  approach  zero.  The  Hamiltonian  is  shown  in  Figure  42.  The  initial 
values  of  the  calculated  costates  are  shown  in  Table  11.  Some  more  details  about  this 
solution  and  required  computer  resources  will  also  be  provided  in  Chapter  III.C.4  when 
comparing  this  solution  with  others. 
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Figure  33. 


Minimum  time  solution  (GPOPS):  3D  plot  of  optimal  rendezvous  trajectory. 
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Figure  34.  Minimum  time  solution  (GPOPS):  2D  plots  of  optimal  rendezvous  trajectory. 
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Figure  35.  Minimum  time  solution  (GPOPS):  control  history  solution. 
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Figure  36.  Minimum  time  solution  (GPOPS):  state  and  costate  histories  for  the  defining 

translational  variables  of  the  chaser  vehicle. 
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Minimum  time  solution  (GPOPS):  state  and  costate  histories  for  the  defining 
angular  parameters  of  the  chaser  vehicle. 
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Minimum  time  solution  (GPOPS):  state  and  costate  histories  for  the  defining 
angular  parameters  of  the  RSO. 
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Figure  38. 
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Figure  39.  Minimum  time  solution  (GPOPS):  time  history  of  the  translational  endpoint 

discrepancies. 


Figure  40.  Minimum  time  solution  (GPOPS):  time  history  of  the  attitude  endpoint 

discrepancies. 

70 


G7 


1 

0.5 
0 

^  -0.5 
-1 
-1.5 


0 


0.5 


y“> 


-0.5 


0  1 


0.5 


0? 


-0.5 


-1 


1 

0 

G7  -1 

-2 

.... 

,  ,  ^  , 

-3 

,  ,  , 

12  3  4 

Time  (s) 


0  12  3  4 

Time  (s) 


Figure  4 1 .  Minimum  time  solution  (GPOPS):  history  of  transversality  conditions. 
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Figure  42.  Minimum  time  solution  (GPOPS):  history  of  the  Hamiltonian. 
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Costate 

Initial  Condition 

Costate 

Initial  Condition 

Costate 

Initial  Condition 

A\ 

-0.0832 

A% 

-0.0300 

A\5 

1.4758 

h 

0.2801 

Ag 

0.0027 

A\6 

1 . 9302 

A?, 

-0 .0172 

-0.1521 

A\n 

0.7304 

A4 

0.2251 

An 

-0.0245 

^-18 

-0.9484 

A5 

-0.4795 

An 

0.5898 

2-19 

-2.1946 

A(, 

0.0426 

^13 

-0.7500 

^  20 

-1.5493 

An 

0.1521 

A\4 

1 .0625 

3 . 4237 

Table  1 1 .  The  initial  values  of  costates  and  tf  as  defined  by  GPOPS  for  minimum  time  solution. 


The  important  thing  about  this  solution  is  that  there  is  no  singular  control  in 
variable  uz  as  there  was  in  the  previous  scenario  described  in  Chapter  II.E.l,  where  the 
translational  controls  were  expressed  in  the  orbital  frame  rather  than  in  the  body  frame. 
However,  the  solution  for  Ty  might  suggest  that  there  may  be  a  singular  arc  in  the 
beginning  of  the  maneuver.  It  should  be  noted  that  in  case  of  a  singular  control,  the 
GPOPS  control  output  may  not  be  feasible  for  onboard  implementation  due  to  it  highly 
oscillating  nature  (Boyarko,  Yakimenko  and  Romano  2009a). 

2.  Minimum  Quadratic  Control  (Minimum  Energy) 

For  the  minimum  energy  solution,  Figures  43  and  44  represent  the  resulting 
trajectory  in  both  three  dimensions  and  planar  view.  Figure  45  shows  a  plot  of  the 
resulting  control  history,  as  well  as  the  associated  switching  function  that  is  used  for 
calculation  of  the  control  based  on  the  running  cost  in  Equation  (61).  Figures  46-48 
show  the  states  and  costates.  The  discrepancies  in  endpoint  conditions  are  shown  to 
reach  zero  at  tf  in  Figures  49  and  50.  The  transversality  conditions  are  shown  in  Figure 
51  and  the  Hamiltonian  is  shown  in  Figure  52.  Table  12  summarizes  the  initial  values  of 
the  costates  that  were  solved  for  by  the  GPOPS  solution.  The  PI  of  the  solution  was  J  = 
0.2445. 
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Figure  44.  Minimum  energy  solution  (GPOPS):  2D  plots  of  optimal  trajectory. 
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Minimum  energy  solution  (GPOPS):  control  history  with  respect  to  the 

optimal  trajectory. 
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Figure  46.  Minimum  energy  solution  (GPOPS):  state  and  costate  histories  for  the 

translational  variables  of  the  chaser  vehicle. 
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Figure  47. 


Figure  48. 
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Minimum  energy  solution  (GPOPS):  state  and  costate  histories  for  the  defining 
angular  parameters  of  the  chaser  vehicle. 
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Minimum  energy  solution  (GPOPS):  state  and  costate  histories  for  the  defining 

angular  parameters  of  the  RSO. 
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Figure  49.  Minimum  energy  solution  (GPOPS):  history  of  the  translational  endpoint 

discrepancies. 


o.i 


-o.i 


<u 

o 

c 

8? 

b 

Q) 

"cc 

OZ 


n  -0.2 

CD 

<  -0.3 


i  i  (  r~ 


"1  i 


0 


Figure  50.  Minimum  energy  solution  (GPOPS):  history  of  the  rotational  endpoint 

discrepancies. 
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Figure  5 1 .  Minimum  energy  solution  (GPOPS):  history  of  the  transversality  conditions. 
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Figure  52.  Minimum  energy  solution  (GPOPS):  history  of  the  Ftamiltonian. 

The  solution  GPOPS  arrives  at  is  feasible  in  that  it  does  not  violate  any  of  the 
constraints  presented  in  the  problem  formulation.  Furthermore,  the  control  history 
follows  the  logic  that  was  derived  from  Minimum  Principle  stated  in  Equation  (23). 


77 


Costate 

Initial  Condition 

Costate 

Initial  Condition 

Costate 

Initial  Condition 

M 

-0.0335 

2-  8 

0.0225 

2-15 

-0.1798 

2-2 

0.0602 

2-9 

0.1005 

2. 16 

0.0457 

h 

-0.0356 

2io 

0.0460 

2-17 

-0.2236 

X4 

-0.1013 

2- 11 

-0.1932 

2-n 

-0.0086 

2-5 

-0.2672 

2- 12 

-0 . 1462 

2*19 

-0.1366 

2-6 

-0 . 0714 

2|3 

0 . 1229 

2  20 

-0.3664 

2-i 

-0.0460 

2l4 

-0.1501 

c 

8 .8943 

Table  12.  The  initial  values  of  costates  and  t  f  as  defined  by  GPOPS  for  minimum  energy 

(quadratic-control)  solution. 


3.  Minimum  Absolute  Control  (Minimum  Fuel) 

To  avoid  using  an  absolute  value  (which  is  a  nondifferentiable  function  at  the 
value  zero),  the  running  cost  in  Equation  (62)  was  substituted  with: 

f0  =  U 4  +  (93) 


with  the  corresponding  modifications  in  the  control  allocation: 


u 


rrri  rri  rri  rri  rri 

lluuulll 

x  max  y  mux  y  max  xmax  y  max  zmax  xmax  j^max  y  max  xmax  ymax  zmax 


(94) 


so  that  now  0  <  u  <  1  (the  superscript  on  the  control,  +  or  -,  denotes  the  direction  in  the 
body  axis).  This  also  results  in  the  following  control  relations  (compare  with  Equation 
(21)): 


fx  (^4  (7l  45  ?3  )(*6  ^7 )  "^  2  (t/l  45  45  45  2  (45  45  +  45  45  )  Ob  ^9) 

fy  =  2  (  4lC  ^2C  +  46  4^  )  («1  -  «7  )  +  (  ?4  2  -  ^1C  2  +  45^  -  <]}2  )  («2  -  «s)  +  2  <?3C  “  Ql  <1*  )  Ob  “  M9  ) 

/r  = 2  ( q"  qCi  -  qci  q* )  Ob  - »7 ) + 2  ( 45c4tf  +  45c45c )  Ob  - u% ) + ( 45“  -  <7C2  -  45C2 + # )  («3  - )  •  (95) 

T=(u4-u]0) 

Ty  =(U5-Un) 

Tz=(u6-un) 


Analyzing  the  structure  of  the  switching  function, 
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H  (M,u) ~XM<  +— (^4/1  +^5/2  +  4;/3)  +  "w 

,_i  /«  7, 


1 


1 


c  /17m4  +  c  4sMs  +  c  4M6 

-'ll  -*22  -*33 


(96) 


and  comparing  it  with  that  of  Equation  (26),  we  can  derive  the  optimal  control  strategy  as 
follows: 


and 


and 


switch,. 

<  —m 

switch,. 

>  —m 

4+3  < 

U 

' — < 

1 

4+3  * 

~!jj 

switch,. 

•6  ^  111 

switch,. 

6  >  m 

4-3 

U 

* — ■ 

VI 

4-3 

>45 

/  =  1,2,3 
i  =  4,5,6,  j  =  i- 3 

i  =  7,8,9 


/  =  10,11,12,  y  =1— 9 


(97) 


Here,  switch,. ,  for  i  =  1,2,3 ,  are  the  same  as  in  Equation  (68).  Hence  in  this  case  we  can 
expect  bang-off-bang  control. 

The  results  are  presented  in  Figures  53-62  and  Table  13.  The  value  of  the  PI  was 
computed  as  J  =  2.4851.  In  Figure  20,  the  controls  ui ,  /  =  1,...,6,  acting  in  the  positive 

direction  along  the  corresponding  body  axis  appear  as  positive  values,  and  uj ,  i  =  7,...,  12 
are  shown  with  the  negative  sign.  Note,  that  according  to  the  optimal  control  structure  of 
Equation  (97)  the  switching  occur  when  the  switches  exceed  a  certain  value,  not  zero  as 
for  example  in  Equation  (69).  For  instance,  for  uz  the  motion  starts  with  some  short 
control  input  in  the  negative  direction,  corresponding  to  u~ ,  and  ends  with  some  short 
control  input  in  the  positive  direction,  corresponding  to  u+  (during  the  remaining  time  the 
control  u.  is  zero).  Also,  the  body  mounted  translational  actuator,  ux,  shows  rapidly 

oscillating  control  values  near  the  end  of  the  maneuver,  which  is  a  property  of  a  singular 
arc.  The  resulting  rapidly  oscillating  control  in  that  region  would  make  the  solution 
incredibly  difficult  to  implement. 
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Figure  53.  Minimum  fuel  solution  (GPOPS):  3D  plot  of  optimal  rendezvous  trajectory. 
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Figure  54.  Minimum  fuel  solution  (GPOPS):  2D  plots  of  optimal  rendezvous  trajectory. 
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Minimum  fuel  solution  (GPOPS):  control  history  with  respect  to  the  optimal 
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Figure  56.  Minimum  fuel  solution  (GPOPS):  state  and  costate  histories  for  the 

translational  variables  of  the  chaser  vehicle. 
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Figure  57. 
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Minimum  fuel  solution  (GPOPS):  state  and  costate  histories  for  the  defining 
angular  parameters  of  the  chaser  vehicle. 
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Minimum  fuel  solution  (GPOPS):  state  and  costate  histories  for  the  RSO. 
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Figure  58. 
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Figure  59.  Minimum  fuel  solution  (GPOPS):  history  of  the  translational  endpoint 

discrepancies. 
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Figure  60.  Minimum  fuel  solution  (GPOPS):  history  of  the  rotational  endpoint 

discrepancies. 
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Figure  61 .  Minimum  fuel  solution  (GPOPS):  history  of  the  transversality  conditions. 
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Figure  62.  Minimum  fuel  solution  (GPOPS):  history  of  the  Ftamiltonian. 
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Costate 

Initial  Condition 

Costate 

Initial  Condition 

Costate 

Initial  Condition 

Ai 

-0.2725 

A& 

-1 . 0009 

A\5 

1 . 6243 

Ai 

0.4399 

a9 

1 .0180 

Ak, 

0.9074 

Ai 

-0.1960 

Aio 

0 .7419 

An 

1.1494 

A4 

0 . 1341 

An 

-0.3327 

An 

-3.2616 

As 

-1.3660 

An 

-0.1168 

Aig 

-A. 2231 

A(, 

-0.8983 

An 

-1.5419 

A-  20 

-2.4355 

Ai 

-0 .7419 

Al4 

2.6216 

f 

8 .1182 

Table  13.  The  initial  values  of  costates  and  tf  as  defined  by  GPOPS  for  minimum  fuel  solution. 


4.  Solution  Comparisons  and  Propagation 
a.  Solution  and  Comparison 

Table  14  summarizes  the  results  for  each  200-node  solution  presented  in 
Chapter  III.C.  Obviously,  the  scenario  with  the  minimum  time  PI  uses  the  most  fuel  to 
complete  the  maneuver.  The  minimum  energy  solution  has  an  order  of  magnitude  lower 
cost,  as  it  spreads  the  control  over  the  entire  maneuver.  It  also  exhibits  a  smoothest 
controls  time  histories,  but  maneuver  takes  the  longest  to  complete.  The  minimum  fuel 
solution  results  in  about  the  same  duration  of  the  maneuver  as  the  minimum  energy 
solution,  but  uses  about  30%  less  fuel.  As  opposed  to  minimum  energy,  it  exhibits  a 
bang-off-bang  control  structure  that  was  expected  with  this  type  of  the  running  cost. 


Minimum  time 

Minimum-Quad- 

Control 

Minimum  fuel 

Time  of  the  maneuver 

3.5086 

8.8943 

8.1182 

Energy  expenditure 

10.4471 

0.2445 

1.1587 

Fuel  expenditure 

20.9419 

3.6941 

2.4863 

Computation  time  (sec) 

33,370 

86,040 

84,261 

Table  14.  Comparison  of  solutions  for  the  three  optimal  control  problems. 


Table  14  also  presents  the  computational  time  it  took  to  converge  to  a 
solution  from  the  very  same  initial  guess.  All  computations  were  performed  on  a 
2.33GHz  Dell  Precision  M90  desktop  computer  with  an  Intel  T7600  processor  and  1Gb 
of  RAM.  As  compared  to  the  results  discussed  in  Chapter  II  (Boyarko  et  al.  2009a),  the 
required  central  processing  unit  (CPU)  time  is  about  three  times  longer.  This  is  due  to  the 
fact  that  200  nodes  were  used  in  the  simulation  as  compared  to  150  nodes  to  that  of 
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Chapter  II.  As  explained  by  Yakimenko  (Yakimenko,  Xu  and  Basset  2008),  the  required 
CPU  time  can  be  drastically  decreased  by  decreasing  the  number  of  nodes,  but  apparently 
the  lesser-node  solution  would  probably  not  catch  short  impulses  of  Figure  47. 

The  ultimate  goal  is  to  allow  trajectory  generation  to  be  carried  in  real 
time  onboard  spacecraft.  At  this  point,  this  cannot  be  realized  using  either  pseudospectral 
methods  (Yakimenko  et  al.  2008). 

b.  Propagating  the  Solution 

The  goal  of  this  section  is  to  check  the  feasibility  of  using  the  optimal 
solution  in  order  to  control  the  system  in  feed-forward  mode.  Consider  the  very  first 
case,  when  the  minimum  time  solution  was  obtained.  Now  that  we  have  the  controls,  let 
us  integrate  the  equation  of  motion  derived  in  Chapter  III.A  with  a  time  step  of  0.0001 
seconds  and  use  the  calculated  controls  as  an  input.  The  resulting  trajectory  is  shown  in 
Figures  63  and  64.  The  time  history  of  the  differences  in  chaser  and  RSO  docking 
position,  velocity,  angular  position  and  angular  rate  are  shown  in  Figure  65.  Table  15  lists 
the  resulting  state  variable  deviations  from  the  desired  endpoints  stated  in  Equations  (55)- 
(58)  at  the  final  time. 

The  propagated  trajectory  appears  to  be  very  close  to  that  of  Figures  32 
and  32,  i.e.,  propagation  of  the  dynamics  using  the  minimal  time  control  output  from 
GPOPS  results  in  a  feasible  trajectory  that  is  similar  to  the  resulting  trajectory  generated 
by  GPOPS.  This  means  that  if  the  solution  could  be  obtained  in  real  time,  with  a  perfect 
model  of  the  real  system,  it  can  be  passed  to  the  control  system  by  feed-forwarding  the 
controls  time  histories. 
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Figure  63.  Propagated  trajectory  using  the  minimum  time  control  history  from  GPOPS. 
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Figure  64.  Propagated  trajectory  using  the  minimum  time  control  history. 
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Figure  65.  Endpoint  conditions  on  chaser  and  RSO  docking  position,  velocity,  angular 

position  and  angular  rate. 


e 

Resulting  Value 

e 

Resulting  Value 

e 

Resulting  Value 

e\ (tf) 

-0.0003 

es  (tf) 

0.0001 

m 

N/A 

ei{tf) 

-0.0021 

e&f) 

-0.0001 

m 

N/A 

e<tf) 

-0.0007 

eio  (tf) 

0.0005 

60 f) 

N/A 

(tf) 

-0.0003 

en  (tf) 

0.0004 

m 

N/A 

es  (tf) 

0.0015 

en(tf) 

-0.0002 

Utf) 

N/A 

ee(tj) 

-0.0014 

en{tf) 

-0.0000 

m 

N/A 

-0.0001 

Ci  (tf) 

N/A 

H(tf) 

N/A 

Table  15.  Value  of  terminal  conditions  at  the  final  time. 
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IV.  TIME-OPTIMAL  REORIENTATION  OF  A  SPACECRAFT 
USING  A  DIRECT  OPTIMIZATION  METHOD  BASED  ON 
QUATERNION  REPRESENTATION  OF  THE  INVERSE 

DYNAMICS 

This  chapter  of  the  manuscript  focuses  on  using  the  equations  of  motion  to  define 
the  controls  for  a  given  attitude  trajectory.  This  is  the  first  step  in  defining  the  necessary 
equations  and  relations  for  the  IDVD  method  for  a  close  approach  that  would  enable 
inspection  or  docking.  In  later  sections,  the  methodology  will  be  applied  to  the 
translational  portion  of  the  problem,  but  first  the  attitude  control  is  considered.  An  initial 
trajectory  is  provided  as  an  initial  guess.  The  trajectory  is  then  perturbed  in  the  solution 
space  by  varying  higher  order  derivatives  on  the  endpoints  until  an  acceptable  solution  is 
found.  It  should  be  noted  that  using  inverse  dynamics  to  optimize  a  rotational  motion  of 
a  satellite  has  been  already  evaluated  by  other  authors  as  well  (Louembet,  Cazaurang, 
Zolghardi,  Charbonnel,  and  Pittet  2007;  and  Mclnnes  1998).  However,  the  Euler  angles, 
suffering  from  well-known  kinematic  singularities  were  used,  and  no  attempt  to  decouple 
the  domains  of  space  and  time  were  made.  Furthermore,  this  research  extends  the 
situations  to  more  realistic  scenarios  where  nonzero  rates  and  accelerations  are  present  as 
in  cases  of  target  tracking  and  docking. 

A.  SPACECRAFT  MODEL  AND  ATTITUDE  TRAJECTORY  OPTIMIZATION 
PROBLEM 

The  first  step  is  to  define  the  dynamics  of  the  system  in  question,  a  reorienting 
spacecraft.  The  rotational  dynamics  of  the  spacecraft  can  be  described  by  Euler’s 
rotational  equations  of  motion.  Written  in  the  body-fixed  principal  axis,  this  results  in  the 
vector  equation  (Greenwood  1987;  Wie  1988): 

Id)  +  (o  x  Ira  =  T ,  (98) 

which  expands  into  the  scalar  equations: 
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(99) 


In  turn,  rotational  kinematics  can  be  described  using  quaternions,  q  =  [c/, ,q2,q3,q4]' .  The 
definition  of  the  quaternion  is  taken  from  Equation  (8)  where  a  is  the  scalar  value  of  the 
rotational  displacement  about  the  eigenaxis  and  p  is  unit  vector  that  describes  the 
direction  of  the  eigenaxis.  The  dynamics  of  the  unit  quaternion  are  described  as 
(Greenwood  1987;  Wie  1988): 
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In  Equations  (99)  and  (100),  I  =  diag([In,I22,I33])  is  the  inertia  matrix  (along  the 
principal  axes),  to  =  [cox ,  co  ,  co.  ]T  is  the  vector  of  angular  velocities,  and  T  =  [Tx ,  Ty ,  7)  ]  ' 
is  the  vector  torques  (bounded  controls). 

The  problem  in  question  is  to  find  the  slew  trajectory  (quaternion  time  history)  for 
a  satellite  subject  to  specific  constraints  that  minimizes  the  time  to  complete  the 
maneuver,  t f  .  This  is  expressed  as  minimizing  the  PI: 


J  =  |  dt , 

o 


(101) 


while  reorienting  a  satellite  from  the  initial  conditions  to0 ,  q0  to  final  conditions  o>  f ,  q  ( 
for  a  system  (99)-(100),  subject  to  constraints  on  controls 


T  <  T  <  T 


(102) 
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Bilimoria  and  Wie  (1993)  formulated  this  problem  for  a  rest-to-rest  maneuver  and 
presented  the  solution  using  an  indirect  method.  They  showed  that,  in  general,  for  the 
case  of  symmetric  body  with  bounds  on  each  torque  component,  it  does  not  result  in  an 
eigenaxis  maneuver.  In  addition  to  that,  the  following  section  presents  even  a  more 
general  solution  obtained  off-line  to  be  used  along  with  that  of  Bilimoria  and  Wie  (1993) 
as  a  reference  one  for  the  proposed  on-line  solution  obtained  using  a  direct  method 
exploiting  inverse  dynamics  of  Equations  (99)-(100). 

Table  1  shows  the  different  test  cases  examined  in  this  chapter.  Test  Case  1  was 
taken  directly  from  Bilimoria  and  Wie  (1993),  while  the  others  were  chosen  to  illustrate 
different  scenarios  of  interest. 


Case 

Normalized  Inertia  Matrix 

Case  1 

l  =  diag([l,l,l]) 

Case  2 

l  =  diag([3,l,2]) 

Table  16. 

Description  of  the  test  cases. 

This  chapter  considers  two  basic  scenarios  assuming  a  0  =90°  and  0=180°  slew 
maneuver  about  the  z-axis  (so  that  q0  =[0,0,0,l]r  and  q;  =  [O,O,sin|0,cosy0]r)  with 
zero  and  nonzero  normalized  body  rates  at  endpoints:  a)  to0  =  to  f .  =  03xl ,  b) 
to0  =  -Q)f  =75-1 3x| ,  and  c)  to0  =  -to^  =  l3xl .  For  the  normalized  states,  the  constraints  (5) 
take  the  form  -l3xl  <  T  <  l3xl . 

B.  SOLVING  THE  PROBLEM  USING  THE  GAUSS  PSEUDOSPECTRAL 
METHOD 

First,  the  problem  formulated  in  Chapter  IV.A  is  addressed  using  the  same 
pseudospectral  direct  method  as  in  Chapters  II  and  III.  The  goal  is  to  have  some  reference 
solutions  for  comparison  to  the  solutions  obtained  through  other  methods,  specifically 
IDVD. 

The  optimal  (Bilimoria  and  Wie  1993)  solution  has  been  matched  using  a 
different  optimization  method  in  previous  work  (Fleming  2004). 
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Figures  66-68  show  the  results  of  applying  GPOPS  to  obtain  minimum  time 
solutions  for  a  180°  slew  of  a  satellite.  Specifically,  Figures  66  and  67  present  time 
histories  of  all  states  and  controls  for  the  solution  that  involves  100  nodes,  which  results 
in  (100-2)x  10=980  variable  parameters,  derived  from  the  98  internal  node  points  for  each 
of  the  7  state  and  3  control  variables.  Figure  68  depicts  the  3D  representation  of  the 
solution  in  inertial  space,  clearly  showing  that  it  is  not  an  eigenaxis  maneuver,  with  an 
inclination  of  the  zb  axis  during  rotation  in  the  xhyh  direction. 

This  solution  compares  with  the  solution  presented  in  Bilimoria  and  Wie  (1993) 
fairly  well.  The  final  calculated  maneuver  time,  tf ,  was  found  to  be  3.243  seconds. 

However,  it  took  almost  two  hours  of  CPU  time  to  obtain  this  solution.  Another 
observation  is  that  because  of  the  nature  of  the  system  described  by  Equation  (99),  the 
optimal  control  has  a  bang-bang  structure  as  shown  in  Figure  67.  That  results  in  the 
maximum  magnitude  of  the  angular  acceleration  at  the  boundary  points  (for  Case  1 
angular  accelerations  are  simply  equal  to  the  corresponding  controls).  It  means  arriving  at 
the  terminal  conditions  with  the  maximum  angular  acceleration.  Also,  if  we  have  to 
update  a  trajectory  while  the  satellite  performs  this  rotation  (to  accommodate  possible 
disturbances  and  unmodeled  dynamics),  it  would  cause  discontinuities  of  angular 
acceleration  (sudden  jumps  in  controls).  Increasing  the  order  of  the  system  to  account  for 
the  boundary  conditions  on  angular  accelerations  will  obviously  cause  a  slight 
degradation  of  the  PI  and  a  further  increase  of  the  required  CPU  time  to  obtain  a  solution. 
Hence,  although  in  this  case  GPOPS  does  produce  a  valid  solution,  as  previously  stated  in 
Chapters  II  and  III,  it  is  not  practical  and  currently  cannot  be  used  for  online 
computations. 
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Figure  66.  Case  1  (GPOPS  solution):  time  histories  of  the  state  variables. 
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Figure  67. 


Case  1  (GPOPS  solution):  time  history  of  the  controlling  torques. 
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Figure  68.  Case  1  (GPOPS  solution):  the  3D  representation  of  the  solution. 

As  pointed  out  in  GPOPS  documentation  (Rao  et  al.  2009),  reducing  the  number 
of  nodes  may  lead  to  a  more  robust  (in  terms  of  computational  time)  result,  therefore  an 
attempt  was  made  to  obtain  a  solution  of  the  same  problem  using  a  lesser  number  of 
nodes.  These  GPOPS  solutions  are  shown  in  Figures  69-71. 

It  turns  out  that  for  a  lesser  number  of  nodes  the  GPOPS  converges  to  different 
solutions.  To  this  end,  Figure  68  shows  time  histories  of  the  angular  velocity  components 
for  the  25-  and  50-node  solution  (involving  230  and  480  varied  parameters,  respectively). 
Obviously,  they  are  different  from  that  of  the  100-node  solution  in  Figure  64.  While  a  25- 
node  solution  is  simply  symmetrical  with  respect  to  the  100-node  solution,  as  can  be  seen 
by  comparing  Figure  66  and  Figure  68  showing  an  inclination  of  the  zh  axis  during 

rotation  in  the  -xbyb  direction,  and  represents  another  equally  optimal  solution  out  of 
possible  four  solutions  (Fleming  2004).  This  is  due  to  the  equivalent  possibilities  of 
positive  or  negative  inclination  in  the  zh  direction  coupled  with  possibilities  of  -xhyh 

and  xbyb  directional  rotation  for  the  180°  maneuver  about  a  principal  axis.  The 
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possibility  of  the  positive  or  negative  inclination  in  the  zb  direction  also  leads  to  two 

equivalent  solutions  for  similar  slews  less  than  180°.  A  50-node  solution,  shown  in 
Figure  69,  appears  not  to  be  valid  (optimal)  at  all. 

As  expected,  decreasing  the  number  of  nodes  leads  to  a  substantial  decrease  in  the 
computational  time,  but  as  shown  above  the  method  could  produce  a  nonvalid  solution. 
Also,  even  if  it  produces  a  valid  solution  the  time  histories  for  control  torques,  Figure  70, 
may  not  be  trackable  by  the  inner-loop  controllers.  Therefore,  the  solutions  obtained  by 
GPOPS  may  not  be  used  in  a  real  time  feed-forward  control  scheme. 


Figure  69.  Case  1  (GPOPS  solution):  comparison  of  time  histories  of  angular  velocity 
components  obtained  for  25  nodes  (top)  and  50  nodes  (bottom). 
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Figure  70.  Case  1  (GPOPS  solution):  comparison  of  time  histories  of  torques,  obtained 

for  25  nodes  (top)  and  50  nodes  (bottom). 


Figure  71.  Case  1  (GPOPS  solution):  the  3D  representation  of  the  25-node  solution. 

For  Case  2,  the  nonsymmetrical  inertia  with  the  bounds  on  individual  control 
torques,  the  solution  is  slightly  different  (Figures  72-73).  The  overall  characteristic  of 
this  and  other  solutions  involving  different  sets  of  the  boundary  conditions  will  be 
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presented  in  Chapter  IV.E,  but  the  general  tendency  is  the  same — it  requires  at  least  a 
hundred  nodes  to  produce  a  valid  and  feasible  off-line  solution.  Yet,  GPOPS  again 
presents  a  good  and  easy-to-use  tool  to  produce  reference  trajectories  that  can  be  used  for 
comparison  with  solutions  obtained  using  other  approaches. 


Figure  72.  Case  2  (GPOPS  solution):  time  histories  of  the  state  variables. 
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Figure  73.  Case  2  (GPOPS  solution):  time  history  of  the  controlling  torques. 
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C.  APPLICATION  OF  INVERSE  DYNAMICS  IN  THE  VIRTUAL  DOMAIN 

APPROACH  WITH  QUATERNION  ATTITUDE  REPRESENTATION 

One  of  the  two  main  ideas  of  Inverse  Dynamics  in  the  Virtual  Domain  (IDVD) 
method  is  exploiting  the  differential  flatness  property  of  the  equations  of  motion 
(Yakimenko  2000;  Cowling,  Yakimenko,  Whidborne  and  Cooke  2007;  Yakimenko  and 
Siegers  2009).  In  the  above  problem,  this  relates  to  the  fact  that  all  the  state  and  control 
variables  can  be  expressed  as  explicit  functions  of  the  output  variable  or  time  derivatives 
of  the  output  variable,  which  in  this  case  is  the  quaternion  itself: 

®=/i(q.q)»  T=/2(q>q>q)-  (103) 

Another  aspect  of  IDVD  involves  handling  computations  in  the  virtual  domain  allowing 
space  and  time  decoupling.  By  doing  so,  a  trajectory  can  be  computed  while  also 
manipulating  the  speed  at  which  that  trajectory  is  followed.  The  following  sections 
present  a  novel  parameterization  for  the  output  variable,  components  of  the  quaternion, 
and  develop  a  step-by-step  computational  scheme. 

1.  Quaternion  Parameterization 

In  order  to  parameterize  the  problem,  the  output  trajectory  is  approximated  using 
some  combination  of  basis  functions.  The  standard  approach  would  be  to  choose  some 
combination  of  polynomials  or  trigonometric  functions  for  the  output  variables 
(Yakimenko  2000;  Cowling  et  al.  2007;  Yakimenko  and  Siegers  2009).  While  this  may 
be  straightforward  when  dealing  with  state  variables  in  translational  space,  it  may 
become  more  challenging  when  dealing  with  expressions  for  orientation. 

While  a  quaternion  may  be  the  preferred  method  to  express  attitude  because  of  the 
lack  of  singularities,  choosing  basis  functions  becomes  more  challenging  because  a 
nonlinear  unit  norm  condition  needs  to  be  preserved  across  the  quaternion  history  (Milam 
2003).  Simply  choosing  arbitrary  polynomials  for  the  1st  three  components  of  the 
quaternion  and  attempting  to  enforce  the  quaternion  constraint  is  insufficient  because  it 
cannot  be  guaranteed  the  constraint  will  be  met  throughout  the  maneuver  (specifically  if 
one  of  the  iterations  on  a  polynomial  provides  a  norm  >  1  in  any  one  of  the  components), 

so  simply  tacking  on  a  penalty  function  for  violating  the  constraint  will  not  work.  For 
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this  reason,  a  specific  polynomial  expression  for  the  quaternion  was  chosen  inspired  by 
the  work  of  Kim,  Kim  and  Shin  (1995).  This  consists  of  expressing  the  quaternion  time 
history  as  an  exponential  function  containing  a  constant  parameter  multiplied  by  a  Bezier 
polynomial  of  a  degree  n: 

q(r)=qoriexP(®A«(r))’  (104) 

i= 1 

where 

A»=X/Ur)  <105> 

j=i 

and 

Pin{r)  =  [n\\ -rrv.  (106) 

A  complete  list  of  quaternion  properties  employed  in  this  dissertation  is  stated  in  the 
Appendix.  Note,  that  in  Equations  ( 1 04)— (1 06)  z  e  [0;  l]  an  abstract  argument  is  used 

instead  of  time.  This  allows  us  to  exploit  certain  attributes  of  the  Bezier  polynomials  and 
define  properties  at  the  beginning  and  endpoints.  For  example,  the  analytic  expressions 

of  the  5th  order  Bezier  Polynomial,  >  are  as  follows: 

/?15(r)  =  r5  -5r4  +  10r3  -10r2  +5r, 

/?25(r)  =  -4r5+15r4-20r3+10r2, 

/?35(r)  =  6r5-15r4  +  10r3,  (107) 

A,5(0  =  -4r5  +5r4, 
p55{z)  =  z\ 

In  this  case, 

q(z-)  =  q0  exp  (to^,  5  (r))  exp  («2^2  5  (r))  exp  (to3^3  5  (r))  exp  (to4^4  5  (r))  exp  (<b5/355  (r)) .( 1 08) 
A  plot  of  the  5th  order  polynomials  is  shown  in  Figure  74. 
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Figure  74.  Value  of  5th  order  Bezier  Polynomials  with  respect  to  the  virtual  arc. 

These  expressions  have  the  favorable  properties: 

A5(0)  =  0  and  Pi  5  (1)  =  1 ,  for  /  =  1,...,5;  (109) 

A5(0)  =  A5(l)  =  5,  As(1)  =  A5(0)  =  0, 

A5(0)  =  A5a)  =  0,  for/- =  2,3,4; 

A5(0)  =  A5(  1)  =  20,  A5(1)  =  A,5(0)  =  0, 

A.5(0)=A,5(i)=-2o,  A,a)=A.5(0)=o,  (in) 

A5(0)  =  A,5(1)  =  0, 

which  fix  the  value  of  the  polynomial  and  its  derivatives  at  the  endpoints  specified  by 
values  of  re  [0;  Tf  j  where  we  set  rf  =  I  that  will  act  as  an  abstract  virtual  time  argument 

to  exploit  the  above  properties  in  Equations  (109)— (1 11).  The  values  of  b)(  are 
coefficients  defined  by  the  following  relation: 

w,  =  Mq^q, ),  for  /  =  1, . . . ,  5  .  (112) 
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Here,  q  are  the  constant  column  vectors  that  act  as  control  points  and  d)(  represents  a 
constant  augmented  angular  velocity  vector  based  on  Equation  (112). 

The  prevailing  idea  is  that  at  q(r  =  0)  =  q0  and  q(r  =  l)  =  q5,  where  q0  and  q5 

can  be  fixed  such  that  q0  =q(t0)  and  q5  =  q ( /  f ) .  This  also  results  in  a  straightforward 

calculation  of  higher  order  derivatives  of  the  quaternion  curve  with  respect  to  the  virtual 
domain  argument  t  .  Results  for  the  first  derivative  were  presented  in  Kim  for  a  3rd  order 
Bezier  polynomial  (Kim  et  al.  1995).  The  results  of  the  2nd  order  derivative  of  a  5th 
order  Bezier  based  quaternion  is  shown  in  Equation  (113): 
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^3  |  ^3 


=  q0  expOS^  (r  ))(<», /?15(r))2  exp(©2^2f5(r))exp(©3^3>5(r))exp(©4^4>5(r))exp(©5^5i5(r)) 

T= 0 

+  q0  exp(«!  A  5  (T))(<aJh5  (r))  exp  ((b7J25  (r))  exp(<»3/?3  5  (r))  exp(<»4/?4  5  (r))  exp(<»5/?5  5  (r)) 

+  q0  exp(o)l/il5  (r))((aj15  (t))  exp(o)2/i25  (r))(o)2/J25  (r))  exp(e)3/i35  (r))  exp(ra4/?45 (r))  exp(o)5/\5(r)) 

+  q0  exp(»,A,5  (*■))(©!  A  ,5  («■))  exp(©2  A2>5  (r))  exp(©3  A  5  (r))(©3  j03>5  (r))  exp(©4>04>5  (r))  exp(©5  A5>5  (r)) 

+  q0  exp(©j  A,s  (*■))(©!  A,5  (r))  exP(«2  A,s  (r))  exp(w3  A  5  (r))  exp(©4A4>5  (r))(©4  A4,5  (T))  exp(ra5  /?5>5  (0) 

+  q0  exp(©jA>5  (r))(oiiy^15  (r))  exp(©2  A2>5  (r))  exp(©3  A  5  (r))  exp(ra4/?45  (r))  exp(©5  /?5  5  (r))(©5  A  5  (r)) 

+  ^(q°  exp(6r A, 5  (r))  exp(w2/?2  5  (r))(<»2/?25  (r))  exp(©3j03>5  (r))  exp(w4/?45  (r))  exp(co5/?55  (r))  J 
+  ^(q°  exp(©j  Aj5  (r))  exp(w2  A2>5  (r))  exp(w3  A  5  (r))(©3A3>5  (r))  exp(co4/?45  (r))  exp(©5  A  5  (r))  J 
+  ^(q°  exp(©jAj5  0))  exp(w2/?25  (r))  exp(©3  /?3  5  (r))  exp(ra4/?4  5  (r))(©4A,5  (0)  exp(w5  A  5  (r))j 
+  ^  (q°  exp(©j  A,s  0))  exp(©2  A2>5  (r))  exp(«3/?35  (r))  exp(ra4/?45  (r))  exp(ra5/?55  (r))(©5  /?5  5  (0))  •  (113) 
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Differentials  of  quaternion  trajectories  based  on  varying  order  Bezier  polynomials 
can  be  analytically  calculated  by  proper  application  of  the  chain  rule.  The  polynomial 
was  expanded  from  a  3rd  to  a  5th  order  in  order  to  fix  the  first  and  second  derivative  of 
the  quaternion  function  at  the  endpoints.  By  applying  Equations  ( 105)— (109),  therefore 

exploiting  the  property  that  certain  terms  /?.  5 ,  /?.  5  and  /)  5  are  equal  to  zero,  derivative 
values  at  the  endpoints  can  be  calculated  by  the  simple  expression: 


dq 

dr  r=0 


5q0<»i, 


dq 

dr  T=l 


5q5d>5 


(114) 


Note,  in  Equation  (114),  the  first  derivative,  how  the  parameter  to,  is  related  to  the 
angular  velocity  at  the  endpoints.  In  similar  fashion, 


d2  q 
~d? 


Ir=0 


=  ^Oq,,®!  +  25q0(b12  +  20q0(O2 , 


d2  q 
dr2 


It=1 


=  -20q4d>4q4  'q5 +20q5w5 +25q5ra52. 


(115) 


2.  Mapping  from  the  Virtual  Arc  to  the  Time  Domain 

Now  that  the  trajectory  is  set  using  a  virtual  domain,  a  mapping  must  be 
employed  to  convert  this  trajectory  into  a  time  dependent  one.  To  do  this,  a  speed 
factor,  X  ,  is  defined  that  maps  the  points  on  the  trajectory  from  the  virtual  domain  to  the 
time  domain,  therefore  defining  the  final  time  of  the  maneuver: 


X  = 


dr 
~dt ’ 


dr 

T 


(116) 


The  Speed  factor  can  be  a  constant  parameter  that  simply  stretches  or  shrinks  time 
evenly  or  it  can  take  the  shape  of  more  complex  function.  For  this  application,  X (r)  is 
restricted  to  a  function  that  contains  a  reduced  number  of  varied  parameters,  developed 
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for  this  dissertation,  that  still  allow  for  the  speeding  up  and  slowing  down  along  a 
trajectory.  The  representation  is  shown  in  Equation  (117). 

A(r)  =  A0  +  At2  +  (1  -  r)2f?  +  (l  -  (1  -  r)2 )  C  +  (1  -  r2)D  (117) 

Obviously,  more  complex  structures  of  A(z)  will  provide  more  flexibility  in  the 
trajectory.  One  main  idea  of  Equation  (117)  is  to  keep  A(z)  positive  for  all  time,  as 
A(z)< 0  would  imply  time  marching  backwards  and  A(z)  =0  signifies  that  time  has 
stopped.  Also,  keeping  the  speed  factor  of  a  form  such  that  an  analytic  integral  to 
Equation  (117)  exists  not  only  provides  computational  efficiency  and  an  accurate 
integration  to  the  minimum  time  PI  but  also  provides  a  continuous  mapping  from  the 
virtual  domain  to  the  time  domain.  Alternatively  stated,  a  continuous  control  history  is 
available,  whose  resolution  does  not  suffer  from  a  limited  number  of  node  points.  This 
attribute  is  later  revisited  in  Chapter  VI.D.2. 

Although  a  speed  vector  of  the  form  Equation  (117)  does  not  allow  matching  the 
optimal  minimum  time  solutions  exactly,  varying  the  parameters  contained  within  A(z) 
(A0,  A,  B,  C,  and  D)  still  allows  variation  of  the  speed  along  the  trajectory  defined  by 
Equation  (104)  to  produce  feasible  and  easy  to  track  suboptimal  solutions. 

3.  Inverting  the  Dynamics 

As  a  result  of  the  mapping  from  virtual  to  time  domain,  the  expression  for  the 
differential  of  q  with  respect  to  time  is: 

(118) 

dt  dr  dt~  dz  dz  dz 

Now,  if  the  trajectory  in  the  virtual  domain,  q(r) ,  is  specified,  along  with  the 
speed  trajectory,  A(z),  the  resulting  trajectory  of  , q(/),  as  well  as  its  higher  order 
derivatives,  can  be  analytically  expressed  and  mapped  to  the  time  domain. 
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Inverting  kinematic  equations  (100)  and  differentiating  the  result  yields  analytical 
expressions  for  the  angular  velocity  and  angular  acceleration: 

<0(0  =  2q_1(0q(0»  (119) 

a(0  =  q_1(02q(0-q(0<o(0- 

The  torque  history  needed  to  follow  such  a  trajectory  is  calculated  by  inverting  Equation 
(99): 

Tx  ( t )  =  ax  ( t )  +  (/33  -  122  )coy  (i t)(oz  (; t ), 

Ty(t)  =  ay(t)  +  (Iu  -I33)ojt (t)ojz(t),  (120) 

Tz  ( t )  =  az  (0  +  (/22  -  /, ,  )®v  (0®,  (0- 


4.  Matching  Endpoint  Conditions 

From  the  preceding  equations,  a  quaternion  history  can  be  developed  based  on  the 
Bezier  polynomial  that  satisfies  predefined  beginning  and  ending  quaternion  values  as 
well  as  setting  the  angular  velocity  and  angular  acceleration  at  the  endpoints.  The  desired 
angular  velocity  and  acceleration  dictate  the  quaternion  derivatives  at  the  endpoints  to  be 
as  follows: 


q(h))  =  q(Ow(h))> 

q(tf)  =  q(tf)a(tf), 

q(t0)  =  q(t0)2a(t0)  +  q(t0)(o(t0), 

qfi/)  =  q(//)2«0/) + q(*/ )<*>(*/•  )• 


(121) 


Based  on  the  properties  of  the  5th  order  Bezier  polynomial,  the  coefficients  of  the 
quaternion  expression  are  then  calculated  by: 


<5,  =IU(a  (122) 

w5  =^t°(^/)5  (123) 


to2 


q01q(to)  +  2Oc)i 
20 


-25U)!2 


(124) 


105 


(125) 


to4 


q4-1  (q(^/)  +  2°qstQi  - 25q5co12 ) qs~’q4 

-20 


Here,  q  is  computed  using  the  second  equation  in  Equation  (119)  and  the 
complementary  q(  parameters  defined  as: 


and 


Finally, 


q0=qOo)  =  q«>L0> 
q5  =q(^/)  =  q(0|r=1, 

fii  =  fioexp^), 
q4  =  q5  exp(tb5)_1. 


q2  =  qj  exp(tb2), 

q3  =  q5(exp(d>4)exp(d>5)) 


d>3  =  ln(q2~'q3). 


(126) 


(127) 


(128) 


The  resulting  benefit  of  this  laborious  fonnulation  is  that  the  attitude  trajectory 
history  of  a  4x1  q  vector,  that  satisfies  the  constraints  of  a  unit  quaternion,  can  be 
specified  by  a  reduced  set  of  parameters.  These  parameters  are  the  initial  and  final 
conditions  on  the  quaternion  itself,  as  well  as  values  of  angular  velocity  and  angular 
acceleration  at  those  endpoints. 


5.  Increasing  the  Polynomial  Order 

More  flexibility  in  the  trajectory  can  be  obtained  by  increasing  the  order  of  the 
Bezier  polynomial  used  in  the  basis  function.  For  the  case  of  a  7th  order  polynomial,  the 
same  structure  as  Equation  (104)  is  employed  but  now  with  n=l .  The  resulting  7th  order 
polynomials  are  shown  in  Figure  75. 
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Figure  75.  Value  of  7th  order  Bezier  Polynomials  with  respect  to  the  virtual  arc. 


This  leads  to  the  introduction  of  q6,q7,  d>6  and  d>7,  which  are  defined  to  be 

consistent  with  previous  definitions  from  Equations  (104)  and  (112).  If  the  values  of 
orientation  and  angular  velocity  at  the  endpoints  are  set,  endpoint  conditions  of  angular 
jerk  as  well  as  angular  acceleration  can  be  used  as  varied  parameters.  Setting  a  specified 
(low)  value  for  the  initial  and  final  jerk  can  be  critical  for  slewing  maneuvers  of  flexible 
spacecraft,  in  order  to  avoid  excitation  of  structural  modes.  To  do  this,  the  third 
derivative  of  the  output  vector  is  calculated  as: 


<73q  o3  d2 q  .  dA  ,  d1  q  dA  ,2  dq  d2A  ,2  dq( d A 

-A  H - —  2  A - A  ^ - - - A  H - —A  +' 


dr 


dr  dr  dr  dr 


dr  dr" 


dr 


dr 


A. 


(129) 


The  new  expressions  for  the  virtual  derivatives  are  also  recalculated,  which  are  analogous 
to  Equations  (114)  and  (115)  but  with  the  addition  of  a  third  derivative  to  accommodate 
the  change  between  5th  and  7th  order  polynomial: 
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d3  q 


dr 3 


=  343q0®i  -882q0tof  -420q0to2  +210q0d)1  +210q0d)3  +  882q0d>1d>2 , 


r=0 


d3  q 


dr 


(130) 


=  343q7d)7  -882q6d)6e)7q61q7  +  210tb5q7  +210q7d)7  +882q7d)7 


T=1 


+  882q0d)1(O2  -420to6q7; 


d2  q 


dr 
d2  q 


=  49q0«j  +  42q0d)2  -  42q0d)1 , 


T= 0 


dr 


=  42q7d)7  +  49q7d)7  -  42co6q7 ; 


r=l 


(131) 


dq 

d r  r=0 

dq 

dr  r=1 


=  7q0«i, 

=  7q7w7. 


(132) 


New  values  for  the  constants  that  fix  the  initial  conditions  of  the  quaternion 
trajectory  can  be  calculated  similar  to  the  5th  order  polynomial,  except  and  extra  step 
needs  to  be  taken  to  accommodate  the  third  derivative  of  q  . 

D.  SOLVING  THE  PROBLEM  USING  IDVD  METHOD 

This  section  presents  the  results  of  using  IDVD  method  with  two  different 
parameterizations  to  obtain  the  minimum  time  solutions  of  the  problem  posed  in  Chapter 
III. A.  The  Matlab  function  fmincon  was  used  to  optimize  the  trajectory  while  varying 
either  the  angular  acceleration  (for  the  5th  order  polynomial)  or  the  angular  acceleration 
and  jerk  (for  the  7th  order  polynomial)  at  the  endpoints  of  the  trajectory.  Together  with 
five  varied  parameters  for  1(f) ,  Equation  (117),  it  yields  1 1  varied  parameters  for  the  5th 
order  polynomial  approximation  and  17 — for  the  7th  order  polynomial  approximation. 
The  constraints  are  that  the  resulting  control  must  obey  Equation  (102)  and  T(t)  cannot 


108 


be  <  0  at  any  instant.  Initial  guesses  for  the  angular  acceleration  and  jerk  are  equal  to 
zero  and  the  initial  guess  of  10'4  is  used  for  all  parameters  contained  in  A(t) . 

1.  IDVD  Solutions  Varying  2nd  Derivative  Results 

Figures  76-79  present  the  results  obtained  when  applying  the  IDVD  with  a 
quaternion  based  on  a  5th  order  Bezier  polynomial  (compare  it  with  the  GPOPS  solution 
presented  in  Figures  66-68),  allowing  variation  of  the  2nd  derivative  at  the  endpoints. 
The  solution  was  evaluated  using  100  nodes  (although,  as  opposed  to  GPOPS  because  the 
solution  is  analytic,  it  would  be  no  difference  running  it  for  a  larger  or  smaller  number  of 
nodes),  resulted  in  slightly  larger  tf  but  took  significantly  less  time  to  compute  at  only 

10.0  seconds. 
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Figure  76.  Case  1  (IDVD  5th  order):  time  history  of  the  quaternion. 
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Figure  77.  Case  1  (IDVD  5th  order):  time  histories  of  the  angular  velocity  and 

accelerations. 


Figure  78.  Case  1  (IDVD  5th  order):  time  history  of  controlling  torques. 

The  major  difference  compared  to  the  GPOPS  solution  is  that  the  controls  do  not 
have  a  bang-bang  nature  anymore.  Again,  it  was  done  intentionally  by  the  choice  of  the 
quaternion  parameterization.  It  occurs  that  when  implemented  in  the  real  time  controller, 
these  controls  will  be  easier  to  track.  Also,  having  different  controlling  torques  at  the 
endpoints  means  having  different  angular  accelerations.  While  for  the  GPOPS  solution, 
the  terminal  angular  accelerations  are  at  the  mercy  of  the  optimization  routine  using 
IDVD  allows  matching  them  with  the  current  accelerations,  which  also  makes  the  control 
algorithm  more  robust.  Figure  78  shows  the  speed  factor,  the  key  element  in  matching 
the  virtual  and  time  domains. 
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Figure  79.  Case  1  (IDVD  5th  order):  mapping  the  virtual  and  time. 

Figure  80  shows  an  outline  of  the  slew  in  inertial  space,  clearly  showing  that  it  is  not  an 
eigenaxis  maneuver. 


Figure  80.  Case  1  (IDVD  5th  order):  principal  axis  outline  of  180°  slew. 
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The  resulting  solution  for  the  same  scenario  using  25  nodes  is  shown  in  Figures  81-84. 


Figure  81. 


Figure  82. 


Figure  83. 
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Case  1  (IDVD  5th  order):  time  history  of  the  quaternion  using  25  nodes. 


Case  1  (IDVD  5th  order):  time  histories  of  the  angular  velocity  and  acceleration 

using  25  nodes. 


Case  1  (IDVD  5th  order):  time  history  of  controlling  torques  using  25  nodes. 
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Figure  84.  Case  1  (IDVD  5th  order):  mapping  the  virtual  and  time  domains  using  25 

nodes. 

The  rate  profdes  for  the  solution  of  25,  50,  and  100  nodes  are  shown  in  Figures  85 
and  86.  Each  case  results  in  a  smooth  control  solution  regardless  of  the  number  of  nodes 
chosen.  This  is  due  to  the  construction  of  the  quaternion  history  as  well  as  the  controls. 
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Figure  85. 


Figure  86. 


Case  1  (IDVD  5th  order):  comparison  of  angular  velocity  time  histories, 
obtained  for  25  nodes  (top)  and  50  nodes  (bottom). 


Case  1  (IDVD  5th  order):  comparison  of  torque  time  histories,  obtained  for  25 
nodes  (top)  and  50  nodes  (middle). 
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As  in  the  case  of  GPOPS  solution  for  Case  2,  the  nonsymmetrical  inertia  matrix 
causes  a  certain  changes  as  compared  to  the  symmetric  matrix  solution.  The  100-node 
IDVD  solution  in  this  case  results  in  4.767s  maneuver  and  requires  about  a  minute  to 
compute  (see  Figures  87-90). 


Figure  87.  Case  2  (IDVD  5th  order):  Time  history  of  the  quaternion. 
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Figure  88.  Case  2  (IDVD  5th  order):  time  histories  of  the  angular  velocity  and 

acceleration. 
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Figure  89.  Case  2  (IDVD  5th  order):  time  history  of  controlling  torques. 


Figure  90.  Case  2  (IDVD  5th  order):  mapping  the  virtual  and  time  domains. 

2.  IDVD  Solutions  Varying  2nd-3rd  Derivative  Results 

For  the  sake  of  comparison,  Figures  91-94  present  the  solution  of  the  same 
problem  using  a  quaternion  based  on  a  7th  order  Bezier  polynomial,  allowing  variation  of 
both  the  3rd  and  4th  derivative  at  the  endpoints.  Although  this  slightly  improves  the  PI,  it 
drastically  increases  the  computational  time.  The  following  section  addresses  this  issue  in 
more  detail. 
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Figure  91. 


Case  1  and  2  (IDVD  7th  order):  time  history  of  the  quaternions  Case  1  (top) 

and  Case  2  (bottom). 
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Figure  92. 
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Case  1  and  2  (IDVD  7th  order):  time  histories  of  the  angular  velocity  and 
acceleration  for  Case  1  (top)  and  Case  2  (bottom). 
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Figure  93. 


Case  1  and  2  (IDVD  7th  order):  time  history  of  controlling  torques  for  Case  1 

(top)  and  Case  2  (bottom). 
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Figure  94.  Case  1  and  2  (IDVD  7th  order):  mapping  the  virtual  and  time  domains  for  Case 

1  (top)  and  Case  2  (bottom)  solutions. 

E.  RESULTS  AND  COMPARISONS 

This  section  presents  a  comparison  of  the  results  obtained  using  the  IDVD  method 
with  those  of  the  GPOPS  method.  It  disregards  the  fact  that  the  results  obtained  with 
GPOPS  for  low  number  of  nodes  are  infeasible  but  rather  concentrates  on  the 
computational  advantages  the  IDVD  approach  has  for  any  number  of  intennediate  points 
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(nodes  in  the  case  of  GPOPS).  To  start  with,  Tables  17  and  18  summarize  the  180°  rest- 
to-rest  slew  maneuver  solutions  for  symmetric  and  asymmetric  inertia  matrix  obtained 
using  GPOPS  and  IDVD,  as  discussed  in  Chapters  IV. C  and  IV.D. 

In  these  tables,  all  results  are  compared  against  the  readily  available  eigenaxis 
maneuver  solution.  First,  it  is  shown  that  the  true  optimal  solution,  obtained  offline 
(Bilimoria  and  Wie  1993),  which  is  not  an  eigenaxis  rotation,  provides  about  8.5%  and 
11%  improvement  of  the  PI,  time  of  maneuver,  for  Case  1  and  Case  2,  respectively. 


As  seen  from  Tables  17  and  18,  the  GPOPS  solution  converging  to  one  of  the 
equally  optimal  solutions  (if  at  all),  assures  about  the  same  gain  in  the  PI  as  the  truly 
optimal  one.  But,  again,  it  takes  significant  computational  time.  Specifically,  the  100- 
node  solution  that  does  converge  and  assures  a  smooth  controls  history  takes  about  an 
hour  to  converge. 


Trajectory  Generation 
Method 

Nodes 

Computational 
Time  (sec) 

cost  (t^) 

% 

improvement 

Eigenaxis 

N/A 

N/A 

3.5449 

O 

'sP 

Optimal  (Bilimoria/Wie) 

N/A 

N/A 

3.2431 

8.51  % 

GPOPS 

25 

15.5 

3.2573 

8.11% 

GPOPS 

50 

200.5 

3.2859 

7.31  % 

GPOPS 

100 

5962.8 

3.2430 

8.52% 

IDVD  5th  order 

25 

3.7 

3.4289 

3.27% 

IDVD  5th  order 

50 

4.8 

3.4373 

3.04% 

IDVD  5th  order 

100 

10.0 

3.4382 

3.01  % 

IDVD  7th  order 

25 

44.0 

3.3756 

4.78% 

IDVD  7th  order 

50 

69.8 

3.3769 

4.74% 

IDVD  7th  order 

100 

121.2 

3.3776 

4.72% 

Table  17.  The  180°  rest-to-rest  slew  maneuver  about  the  z-axis,  symmetric  inertia  (Case  1). 
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Trajectory  Generation 
Method 

Nodes 

Computational 
Time  (sec) 

cost  (ff) 

% 

improvement 

Eigenaxis 

N/A 

~0 

5.0133 

~0% 

GPOPS 

25 

22.2 

4.4609 

1 1 .02% 

GPOPS 

50 

520.1 

4.4499 

1 1 .24% 

GPOPS 

100 

1893.9 

4.4528 

11.18% 

IDVD  5th  order 

25 

4.8 

4.7857 

4.54% 

IDVD  5th  order 

50 

7.6 

4.7861 

4.53% 

IDVD  5th  order 

100 

9.5 

4.7864 

4.53% 

IDVD  7th  order 

25 

24.3 

4.6768 

6.71  % 

IDVD  7th  order 

50 

41.2 

4.6846 

6.56% 

IDVD  7th  order 

100 

80.4 

4.6869 

6.51  % 

Table  1 8.  The  1 80°  rest-to-rest  slew  maneuver  about  the  z-axis,  asymmetric  inertia  (Case  2). 


As  discussed  in  the  previous  section,  the  IDVD  solution  has  a  much  more  robust 
performance,  allowing  computing  the  same  type  of  maneuvers  just  in  a  few  seconds  as 
opposed  to  hours.  If  an  executable  optimization  library  is  implemented,  the  IDVD 
method  produces  solutions  in  fractions  of  a  second  (Yakimenko  and  Siegers  2009).  Of 
course,  some  of  the  optimality  (PI  value)  is  sacrificed.  On  the  positive  side,  the  solution 
is  always  feasible  and  smooth  for  any  number  of  computational  points,  and  can  be 
brought  closer  to  the  GPOPS  solution  (in  terms  of  the  value  of  a  PI)  by  increasing  the 
number  of  varied  parameters  (the  order  of  the  quaternion  approximation  polynomial). 
Furthermore,  since  the  resulting  IDVD  solution  is  analytic  in  nature,  increasing  the  node 
points,  after  a  solution  is  obtained,  is  a  trivial  evaluation  of  an  analytic  expression. 

As  discussed  in  Chapter  IV.B,  this  solution  features  a  bang-bang  control,  i.e.,  does 
not  account  for  controllers’  dynamics,  and  therefore  can  still  not  be  used  onboard  as  is. 
On  the  other  hand,  the  always-feasible  and  ready-to-go  IDVD  solution  (employing  as  low 
as  say  25  computational  points)  can  be  produced  much  faster,  but  surrenders  up  to  2A  of 
its  gain  as  compared  to  that  of  the  GPOPS  solution  (about  Vi  for  the  7th  order 
approximation). 

Tables  19  and  20  present  similar  data  for  the  90°  rest-to-rest  slew  maneuver. 
While  GPOPS  provides  about  3%  gain  compared  to  a  simple  eigenaxis  slew  solution,  the 
IDVD  method  has  almost  no  advantage  or  may  be  even  worse  if  using  a  5th  order 

quaternion  approximation.  This  is  because,  at  some  point,  the  gains  made  by  having  the 
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ability  to  calculate  a  bang-bang  solution  (with  discontinuous  controls)  outweighs  the 
benefits  of  having  the  ability  to  exploit  the  gains  made  by  the  effect  of  slewing  off  the 
eigenaxis.  All  major  conclusions,  however,  remain  the  same. 


Trajectory  Generation 
Method 

Nodes 

Computational 
Time  (sec) 

cost  (ff) 

% 

improvement 

Eigenaxis 

N/A 

N/A 

2.5066 

~0% 

GPOPS 

25 

20.7 

2.4336 

2.91  % 

GPOPS 

50 

120.0 

2.4332 

2.93% 

GPOPS 

100 

236.2* 

2.4296 

3.07% 

IDVD  5th  order 

25 

5.8 

2.5654 

-2.34% 

IDVD  5th  order 

50 

7.2 

2.5666 

-2.39% 

IDVD  5th  order 

100 

9.4 

2.5671 

-2.41  % 

IDVD  7th  order 

25 

39.7 

2.5043 

0.09% 

IDVD  7th  order 

50 

58.9 

2.5058 

0.03% 

IDVD  7th  order 

100 

77.3 

2.5058 

0.03% 

Table  19.  The  90°  rest-to-rest  slew  maneuver  about  the  z-axis,  symmetric  inertia  (Case  1). 


Tables  21  and  22  compare  GPOPS  and  IDVD  solutions  for  such  case,  when  other 
sets  of  nonzero  boundary  conditions  were  explored  as  well,  and  proved  to  maintain  the 
same  trends.  In  this  table,  all  results  are  compared  against  a  valid  100-node  GPOPS 
solution.  As  seen,  the  GPOPS  solutions  with  a  lesser  number  of  nodes  produce 
somewhat  infeasible  solutions,  meaning  that  they  cannot  be  implemented  in  the  control 
scheme  explicitly.  The  IDVD  solutions  may  yield  to  GPOPS  as  much  as  about  4%  with 
respect  to  the  PI,  but  again  are  produced  much  faster.  Furthermore,  as  shown  in  Figures 
95  and  96,  the  90°  and  180  0  rest-to-rest  slew  maneuvers  with  zero  boundary  rates  feature 
multiple  equally  optimal  solutions,  so  that  both  GPOPS  and  IDVD  solutions  converge  to 
different  solutions  when  changing  the  number  of  nodes  (GPOPS)  /  computational  points 
(IDVD).  In  contrast,  for  the  case  of  nonzero  boundary  rates,  they  all  converge  to  the 
same  solution, as  shown  in  Figure  97. 
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Trajectory  Generation 
Method 

Nodes 

Computational 
Time  (sec) 

cost  ( tf ) 

% 

improvement 

Eigenaxis 

N/A 

~0 

3.5449 

~0% 

GPOPS 

25 

17.4 

3.4450 

2.82% 

GPOPS 

50 

168.0 

3.4408 

2.94% 

GPOPS 

100 

1432.1 

3.4430 

2.87% 

IDVD  5th  order 

25 

5.9 

3.6277 

-2.33% 

IDVD  5th  order 

50 

9.8 

3.6307 

-2.42% 

IDVD  5th  order 

100 

17.3 

3.6316 

-2.44% 

IDVD  7th  order 

25 

41.5 

3.5373 

0.21  % 

IDVD  7th  order 

50 

80.4 

3.5389 

0.17% 

IDVD  7th  order 

100 

158.8 

3.5394 

0.16% 

Table  20.  The  90°  rest-to-rest  slew  maneuver  about  the  z-axis,  asymmetric  inertia  (Case  2). 


It  should  be  noted  that  in  practice,  the  direct  methods  would  likely  be  used  in 
situations  where  the  end-conditions  (angular  rates,  accelerations)  of  the  slew  are  specified 
and  not  equal  to  zero  (to  meet  mission  requirements  of  matching  attitude  rates  of  a 
tumbling  vehicle,  for  example).  For  this  case,  no  simple  eigenaxis  slew  solution  exists 
and  therefore  any  solution  produced  on-line  would  be  good.  The  next  chapter  of  this 
dissertation  exploits  the  IDVD  and  coupling  with  translational  maneuvers. 


Trajectory  Generation 
Method 

Nodes 

Computational 
Time  (sec) 

cost  (fi) 

%  within 
GPOPS  100 
nodes 

GPOPS 

25 

48.3 

2.4028 

0.07% 

GPOPS 

50 

333.9 

2.4016 

0.02% 

GPOPS 

100 

1182.6 

2.4011 

0.00% 

IDVD  5th  order 

25 

5.9 

2.5437 

5.94% 

IDVD  5th  order 

50 

5.6 

2.5450 

5.99% 

IDVD  5th  order 

100 

6.1 

2.5451 

6.00% 

IDVD  7th  order 

25 

27.9 

2.4885 

3.64% 

IDVD  7th  order 

50 

41.5 

2.4895 

3.68% 

IDVD  7th  order 

100 

90.6 

2.4896 

3.69% 

Table  2 1 .  The  90°  maneuver  for  symmetric  inertia  (Case  1)  and  nonzero  boundary  rates, 

e>0  =  ~(of  =  [0.1  0.1  0.l]T . 
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Trajectory  Generation 
Method 

Nodes 

Computational 
Time  (sec) 

cost  (tf) 

%  within 
GPOPS  100 
nodes 

GPOPS 

25 

11.7 

2.9054 

0.03% 

GPOPS 

50 

111.3 

2.9047 

0.01  % 

GPOPS* 

100 

164.3 

2.9045 

0.00% 

IDVD  5th  order 

25 

6.3 

3.0060 

3.49% 

IDVD  5th  order 

50 

5.5 

3.0064 

3.51  % 

IDVD  5th  order 

100 

8.3 

3.0066 

3.52% 

IDVD  7th  order 

25 

31.4 

2.9559 

1 .77% 

IDVD  7th  order 

50 

60.1 

2.9568 

1 .80% 

IDVD  7th  order 

100 

106.5 

2.9571 

1 .81  % 

Table  22. 


The  90°  maneuver  for  symmetric  inertia  (Case  1)  and  nonzero  boundary  rates, 

CO0  =  — CDy  =  —  [l  1  l]T. 
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Figure  95.  Case  1 :  1 80°  rest-to-rest  slew  profile  as  projected  onto  the  x-y  inertial  plane. 
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Figure  96.  Case  1 :  90°  rest-to-rest  slew  profile  as  projected  onto  the  x-y  inertial  plane. 
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Figure  97.  Case  1 :  90°  slew  profde,  co0  =  -to  f  =  [0. 1  0.1  0.  l]T ,  as  projected  onto  the  x-y  inertial  plane. 
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V.  RAPID  ATTITUDE  AND  TRANSLATIONAL  TRAJECTORY 

GENERATION 


Inverse  Dynamics  in  the  Virtual  Domain  exploits  the  ability  to  take  a  polynomial 
of  sufficiently  high  order  and  match  desired  values,  and  potentially  derivative  values,  at 
its  endpoints  to  generate  a  trajectory.  This  is  done  by  representing  each  position  variable 
of  the  trajectory  of  the  close  approach  maneuver  in  the  orbital  frame  by  a  summation  of 
polynomials  and  trigonometric  basis  functions.  While  some  endpoint  conditions  are 
specified  based  on  the  problem  statement,  higher  order  derivative  of  the  endpoints  can 
function  as  parameters  that  may  be  varied  to  optimize  the  trajectory  based  on  a  chosen  PI. 
The  remainder  of  this  chapter  explains  the  process  in  detail  and  applies  the  IDVD  method 
to  the  benchmark  problem  stated  in  Chapter  III. 

A.  INVERSE  DYNAMICS  IN  THE  VIRTUAL  DOMAIN  FORMULATION 

First,  a  basis  function  for  the  translational  trajectory  is  defined  by  Equation  (133). 


r7i 


n 


p(r)  =  a0  +alr  +  a2r “  +  a2 r  +b0  sin(— r)  +  c0  cos(—  r)  +  bl  sin^rj  +  q  cos(;zt)  .  (133) 


This  leads  to  3  separate  equations  (and  separate  set  of  coefficients)  for  the  jc,  y,  and  z 
position  of  the  chaser  vehicle  each  based  on  the  form  of  the  basis  function  described  by 
Equation  (133).  It  should  be  reinforced  that  the  trajectory  described  by  Equation  (133) 
with  higher  order  derivatives  shown  by  Equations  (134)— (135),  are  constructed  with 

respect  to  a  virtual  argument,  re^Gfz^J,  and  not  time.  This  specifies  the  trajectory  in 

the  spatial  domain,  but  lets  the  speed  the  trajectory  is  traversed  to  be  varied  as  well. 


P  = 


dp 

dv 


f  7T  77  77  77  ^ 

o,  +2a2r  +  3a2r  +b0  —cos (—t)-cq  ysin .(— r) 


+bx7t  cos(tzt)  -cxn  sin(^r) 


J 


dp 

~df 


^  7T1  .  71  TC1  77  ^ 

a2  +6  a3r-b0  —  sin(—  r)  +  c0  —  cos(—  r) 
.2 


-bx7T  sin(^r)  -cpr  cos(/zr) 


J 


(134) 


(135) 
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The  plots  of  the  basis  functions  considered  are  shown  in  Figure  98: 
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Figure  98.  Basis  functions  considered  for  translational  trajectory  generation. 


1.  Mapping  the  Virtual  Arc  and  Specifying  Endpoint  Values 

The  same  structure  of  the  virtual  time  argument  is  used  as  in  Chapter  IV. B. 2, 
Equation  (117).  This  allows  for  5  varied  parameters  that  determine  the  history  of  the 

speed  factor,  ^ ,  over  the  duration  of  the  maneuver.  The  endpoint  values  of  the  basis 
functions  are  simply  set  to  the  beginning  and  ending  values  of  the  maneuver  in  the  spatial 
domain.  The  derivative  values  at  the  endpoints  and  their  derivatives  in  the  virtual  domain 
are  calculated  based  on  the  time  derivatives  at  the  endpoints  and  speed  factor  shown  in 
Equations  (136)  and  (137). 


dp  jit  f) 
dr 


=  Pi(tf)/  Mtf)  , 


dPj(t0) 

dr 


for  i  =  x,y,z 


Pido)  /  Wo) 


(136) 
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and 


d2Pi(tf)  dZ(tf ) 

d  PM o)  \  ■  u  \  dMt0) 


(137) 


dr2 


■  =  (/W-A(*o) 


dt 


)U(t0y- 


In  order  to  specify  the  2nd  derivatives  at  the  endpoints,  a  basis  function  with  a  reduction 
in  the  number  of  coefficients  is  used  (Yakimenko  2000;  Yakimenko  and  Siegers  2009): 


,7t 


p(r)  =  a0  +  axz  +  a2r~  +  c0  cos(—  z)  +  b\  sin(^r)  +  cx  cos(/zr)  . 


(138) 


The  resulting  coefficients  are  then  be  derived  according  to  the  algebraic  relation  in 
Equation  (139). 


“l 

0 

0 

0 

0 

1 

a0 

pit  o) 

1 

1 

1 

1 

0 

-1 

a\ 

p(tf) 

0 

1 

0 

nil 

71 

0 

a2 

p\t o) 

0 

1 

2 

0 

-n 

0 

bo 

p'(tf) 

0 

0 

2 

0 

0 

-n2 

bx 

p"(t  o) 

0 

0 

2 

-n2  /  4 

0 

n2 

_clJ 

p"(tf)_ 

For  higher  order  basis  functions  that  are  capable  of  specifying  jerk  at  the 
endpoints,  the  basis  function  is  specified  as  Equation  (133)  and  the  coefficients  are 
obtained  by  solving  the  algebraic  relation  in  Equation  (140). 
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2.  Implementation  of  Quaternion  Based  Attitude  Representation 

The  attitude  implementation  of  IDVD  for  the  close  rendezvous  problem  described 
in  Chapter  III  is  done  similar  to  Chapter  IV.B.  The  proper  Bezier  polynomial  is  chosen 
such  that  the  desired  boundary  conditions  can  be  met  as  well  as  allowing  variation  of  the 
next  highest  derivative.  One  major  difference  is  that  the  attitude  trajectory  is  formulated 
in  the  orbital  frame  instead  of  the  inertial.  The  angular  rates  with  respect  to  the  inertial 
frame  are  then  computed,  based  on  Equation  (141),  and  used  to  detennine  the  resulting 
control  torques  in  the  body  frame. 

COx  —  OJ  .  +  2 (£/,£/,  +  t/2<74  )f- 

0\  =  °(Oy  +  2(q2q2  -qxqA) Q  (141) 

co.  —  co.  +  2 (q4  —  ql  —  q2  +  q2  )Q 

3.  Alternative  Euler  Angle  Formulation  for  Attitude  Representation 

Considering  an  Euler  angle  representation  of  a  2-1-3  with  respect  to  the  orbital 

reference  frame  rotation  can  be  expressed  by  the  angular  displacements  an<^  ,(// 

pertaining  to  angular  displacements  about  the  y,  x  and  then  z  body  axis  respectively  (Sidi 
1997).  The  angular  rates  with  respect  to  the  orbit  frame  can  then  be  expressed  as: 

°cox  =  (j)  cos  {y/)  +  6  sin(^)  cos(^) 

°  co y  =  -(j)  sin(^)  +  0  cos(^)  cos(^)  (1 42) 

°coz  =\j/  -0  sin(^) 

The  angular  rates  with  respect  to  the  inertial  frame  can  then  be  expressed  as: 

cox  =^cos(^)  +  ^sin(^)cos(^)  +  2(^1^3  +^2^4)Q 

coy  =  -^sin(^-)  +  ^cos(^)cos(^)  +  2(^2^3 -qxqA)Cl  (143) 

co.  =  \j/  -  6 sin(^)  +  2(^f42  -  q ,2  -q22  +q2)Cl 

The  resulting  angular  accelerations,  with  respect  to  the  inertial  frame  expressed  in  the 
body  frame  becomes: 
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ax  =  (f>  cos  <y )  -  (j)\/ 7  sin(^)  +  9  sin(^)  cos(^)  +  9\j/  cos  {(/))  cos(^) 

-  d(j)  sin(^)  sin(^)  +  2(q3qx  +  qxq 3  +  q4q2  +  q2q4)coo 

a  =  -0sin(i//)-<f>i//cos(i//)  +  9cos(i//)cos(<f>)-9i//cos(<t>)sm(i//)  (144) 

-  6(j)  sin(^)  cos(^)  +  2(q2q3  +  q3q2  -  q4qx  -  qxq4  )a>0 

oc,  =y/~  d<j)  cos {(/))  -  9  sin(^)  +  2 (q4q4  -qxqx-  q2q2  +  q2q2 )a>0 

This  equates  to  the  required  inertial  torque  expressed  in  the  body  frame  dictated  by: 

Tx  =Al«v+(/33  -/2 2K®v 

Ty  =  /22«  +(/|,  -4)®^  (145) 

Tz  —  +(^22  ^1 1 

This  formulation  allows  complete  description  of  the  attitude  trajectory  in  the  orbital  and 
inertial  frames  and  the  associated  rates  and  inertial  torques  by  specifying  initial  and  final 
conditions  on  9,  t//,  and  (f)  as  well  as  their  higher  order  derivatives. 


B.  SOLVING  THE  IDVD  PROBLEM 

For  the  close  rendezvous  problem  with  the  IDVD  formulation,  the  Sparse 
NOlinear  OPTimizer  (SNOPT)  (Gill,  Murray,  and  Saunders  1996)  was  used  to  solve  for 
the  resulting  trajectory.  For  this  implementation,  the  objective  function  was  the  PI  based 
on  the  specific  type  of  cost  function.  The  problem  constraints  were  programmed  as: 


control  constraint:  -u  .  <  u  <  um!iv 

min  iiiax 

path  constraint:  ( x 2  +  y2  +z2)-r2  >  0 

time  constraints:  A(t)  >0  Vr 


(146) 


with  the  limits  on  u  and  value  of  r  set  to  those  used  in  Chapter  III.  It  is  reinforced  that, 
with  the  IDVD  method,  the  dynamic  constraints  (equations  of  motion)  of  the  problem  and 
desired  endpoint  conditions  are  always  satisfied  exactly  at  every  iteration.  If  the  problem 
is  not  stated  such  that  a  solution  does  not  exist,  judicious  (or  sometimes  common  sense) 
selection  of  the  initial  guess  on  the  varied  parameters,  conservative  guesses  on  the 
endpoint  derivatives  and  speed  factor  coefficients,  will  lead  to  always  having  a  feasible 
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solution  to  implement.  This  is  not  the  case  with  pseudospectral  methods,  since  at  any 
given  iteration,  the  kinematic  and  dynamic  constraints  in  the  equations  of  motion  are  not 
guaranteed  to  be  satisfied  and  may  result  in  an  infeasible  solution. 

C.  COMPUTATIONAL  RESULTS 

Both  IDVD  methods  of  matching  endpoint  translational  and  angular  position  and 
velocities  are  used  to  generate  solutions.  The  1st  set  of  endpoint  parameters  is  considered 
to  be  the  position  and  the  2nd  set  of  endpoint  parameters  is  the  velocity.  The  IDVD  (3rd) 
method  refers  to  varying  the  set  of  3rd  parameters,  accelerations,  at  the  endpoints  while 
the  IDVD  (3rd-4th)  method  refers  to  varying  both  the  third  and  fourth  parameters, 
accelerations  and  jerk. 

1.  Minimum  Energy  Cost  Varying  3rd  Parameter  Set  Results 

The  results  for  using  the  IDVD  (3rd)  method  for  a  minimum  energy  cost, 
Equation  (15),  are  shown  next.  They  are  based  on  the  translational  basis  Equation  (139) 
and  the  quaternion  attitude  formulation  employing  the  5th  order  Bezier  polynomial.  The 
result  is  suboptimal  when  compared  to  the  infinite  dimensional  optimal  control  problem 
solution,  but  is  optimal  based  on  the  additional  constraint  of  only  using  the  specified 
polynomial  set  of  basis  functions.  This  reduces  the  number  of  varied  parameters  to  17 
(the  third  derivative  of  each  state  plus  the  coefficients  of  the  speed  factor),  resulting  in  a 
computational  time  of  12.2818  seconds.  The  final  PI  based  on  the  minimum  quadratic- 
control  cost  is  J  =  0.2461.  Figure  99  shows  the  3D  representation  of  the  trajectory  in 
Figure  100  shows  the  planar  projection  of  the  same  trajectory.  The  control  solution  is 
shown  in  Figure  101,  which  is  smooth  over  the  interval,  consistent  with  the  optimal 
solution  found  in  Chapter  III.C.2.  Figures  102  and  103  show  the  time  history  of  the 
endpoint  discrepancies  and  verify  that  the  endpoint  conditions  were  met  at  the  final  time. 
The  solution  using  this  specific  formulation  and  cost  is  highly  attractive  due  to  it’s 
inherent  low  fuel  cost  and  the  attribute  that  the  solution  does  not  demand  that  the 
actuators  be  saturated  at  any  given  instant,  unlike  the  bang-bang  and  bang-off-bang 
nature  of  minimum  time  and  minimum  fuel  cost  functions  respectively.  Specifically,  this 
noncontrol-saturated,  smooth  trajectory  is  desirable  as  the  spacecraft  is  in  the  final  stages 
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of  matching  docking  points.  The  speed  factor,  a  representative  measure  of  how  fast  the 
chaser  is  moving  along  the  spatial  trajectory,  is  shown  in  Figure  104. 


Figure  99.  Minimum  energy  solution  (IDVD  3rd):  3D  optimal  rendezvous  trajectory. 
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Figure  101.  Minimum  energy  solution  (ID YD  3rd):  control  history. 


Figure  102.  Minimum  energy  solution  (IDVD  3rd):  history  of  the  translational  endpoint 

discrepancies. 
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Figure  103.  Minimum  energy  solution  (IDVD  3rd):  time  history  of  the  rotational  endpoint 

discrepancies. 
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Figure  104.  Minimum  energy  solution  (IDVD  3rd):  time  history  of  the  speed  factor. 
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2.  Minimum  Time  Cost  Varying  3rd  Parameter  Set  Results 

The  final  PI  based  on  the  minimum  time  cost  is  J  =  3.734  seconds.  Figure  105 
shows  the  3D  representation  of  the  trajectory  in  Figure  106  shows  the  planar  projection  of 
the  same  trajectory.  The  control  solution  is  shown  in  Figure  107,  which  although  is 
smooth  by  construction,  approaches  the  bang-bang  nature  of  an  optimal  solution 
computed  with  respect  to  a  minimum  time  cost  function  presented  in  Chapter  III.C.l. 
Figures  108  and  109  show  the  time  history  of  the  endpoint  discrepancies  and  verify  that 
the  endpoint  conditions  were  met  at  the  final  time.  Figure  110  shows  the  speed  factor, 
which  revels  that  the  chaser  is  speeding  up  along  the  spatial  trajectory  and  then  slowing 
down  as  it  approaches  the  final  desired  states.  This  type  of  qualitative  behavior  is 
inherent  to  bang-bang,  minimum  time  solutions  to  general  maneuvers  which  apply  max 
control  at  the  beginning  to  increase  speed  of  the  maneuver,  then  max  control  at  the  end  to 
decrease  speed  and  arrive  at  the  desired  final  states. 


Figure  105.  Minimum  time  solution  (IDVD  3rd):  the  3D  optimal  rendezvous  trajectory. 
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1 


Figure  106. 


Minimum  time  solution  (ID YD  3rd):  2D  planar  projection  of  the  trajectory. 


Figure  107.  Minimum  time  solution  (ID VD  3rd):  control  history. 
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Figure  108.  Minimum  time  solution  (IDVD  3rd):  history  of  the  translational  endpoint 

discrepancies. 
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Figure  109.  Minimum  time  solution  (IDVD  3rd):  history  of  the  rotational  endpoint 

discrepancies. 
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Figure  110.  Minimum  time  solution  (IDVD  3rd):  time  history  of  the  speed  factor. 

3.  Minimum  Fuel  Cost  Varying  3rd  Parameter  Set  Results 

The  final  PI  based  on  the  minimum  fuel  cost  is  J  =  2.83 12.  Figure  111  shows  the 
3D  representation  of  the  trajectory  in  Figure  112  shows  the  planar  projection  of  the  same 
trajectory.  The  control  solution  is  shown  in  Figure  113.  Upon  close  inspection,  the 
control  solution  for  each  component  has  more  action  at  the  beginning  and  ends  of  the 
maneuver,  exhibiting  qualities  of  a  bang-off-bang  structure  consistent  with  the  cost 
function  presented  in  Chapter  III.C.l.  Figures  114  and  115  show  the  time  history  of  the 
endpoint  discrepancies  and  verify  that  the  endpoint  conditions  were  met  at  the  final  time. 
Figure  116  shows  the  speed  factor,  which  is  similar  to  the  minimum  energy  solution 
presented  in  the  previous  section. 
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Figure  111.  Minimum  fuel  solution  (IDVD  3rd):  the  3D  optimal  rendezvous  trajectory. 
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Figure  1 12.  Minimum  fuel  solution  (IDVD  3rd):  2D  planar  projection  of  the  trajectory. 


142 


0.5 


Figure  1 13.  Minimum  fuel  solution  (IDVD  3rd):  control  history. 


Figure  1 14.  Minimum  fuel  solution  (IDVD  3rd):  history  of  the  translational  endpoint 

discrepancies. 
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Figure  115.  Minimum  fuel  solution  (IDVD  3rd):  history  of  the  rotational  endpoint 

discrepancies. 
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Figure  116.  Minimum  time  solution  (IDVD  3rd):  time  history  of  the  speed  factor. 
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4.  Minimum  Energy  Varying  3rd-4th  Parameter  Set  Results 

The  problem  stated  in  Chapter  V.B.l  is  also  solved  using  the  IDVD  (3rd-4th) 
method  based  on  the  minimum  energy  cost  function  from  Equation  (15).  The  results  are 
extremely  similar  to  the  solution,  as  shown  in  the  3-D  trajectory  of  Figure  117,  to  the 
same  problem  formulation  based  on  the  IDVD  (3rd)  method,  except  the  increased 
flexibility  in  the  polynomial  due  to  the  extra  basis  functions  allows  for  smoother  controls 
at  the  endpoints,  as  shown  in  Figure  118  and  later  in  Figure  123,  when  compared  with 
IDVD  (3rd).  The  solution  took  623  seconds  to  compute  and  had  a  PI  of  J=  0.2458,  only  a 
0.12%  reduction  in  PI.  It  should  be  noted  that  using  this  formulation,  you  can  specify 
constraints  on  the  jerk  profiles  while  still  using  the  acceleration  profile  as  your  control 
vector.  This  attribute  is  based  on  the  polynomial  formulation  of  the  trajectory  and  cannot 
be  implemented  using  pseudospectral  techniques.  To  set  constraints  on  jerk  using 
pseudospectral  techniques,  the  control  vector  would  have  to  be  based  on  jerk  and  not 
acceleration. 


145 


CO 

3 


Figure  118.  Minimum  energy  solution  (ID VD  3rd-4th):  control  history. 


5.  Minimum  Time  Varying  3rd-4th  Parameter  Set  Results 

The  minimum  time  solution  using  the  IDVD  (3rd-4th)  method  is  shown  in  Figure 
119.  It  has  a  PI  of  J  =  .3.6582  and  took  86.6  seconds  to  compute.  The  3D  trajectory 
looks  similar  to  that  obtained  by  IDVD  (3rd),  but  the  extra  flexibility  in  the  polynomial 
provided  by  an  increased  set  of  basis  function  allows  for  shaper  transitions  in  the 
controls,  shown  in  Figure  120,  mimicking  that  of  a  bang-bang  controller.  Although  the 
computational  cost  rose  to  86.6  seconds,  it  is  still  well  below  the  8+  hour  computational 
cost  of  GPOPS  for  a  200  node  representation. 
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Figure  119.  Minimum  time  solution  (IDVD  3rd-4th):  3D  optimal  rendezvous  trajectory. 
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Figure  120.  Minimum  time  solution  (IDVD  3rd-4th):  control  history. 
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6.  Minimum  Fuel  Varying  3rd-4th  Parameter  Set  Results 

The  minimum  fuel  solution  obtained  using  the  IDVD  (3rd-4th)  method  is  shown 
in  Figure  121.  Again,  the  trajectory  is  very  similar  to  the  solution  calculated  from  the 
IDVD  (3rd)  except  the  increased  flexibility  allows  for  shaper  transitions  in  control  at  the 
endpoints.  The  bang-off-bang  type  of  control  is  evident  from  examining  the  control 
histories  in  Figure  122.  The  resulting  PI  is  J  =  2.7552,  taking  405.5  seconds  to  compute. 
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Figure  122.  Minimum  fuel  solution  (IDVD  3rd-4th):  control  history. 

7.  Propagated  IDVD  Solutions 

As  with  the  solutions  obtained  by  the  pseudospectral  method  in  Chapters  II  and 
III,  the  solutions  for  the  controls  obtained  with  the  IDVD  method  were  propagated  using 
the  same  mechanism  as  the  pseudospectral  method  stated  in  Chapter  III.C.4.  The 
resulting  discrepancies  in  the  endpoints  resulting  from  all  three  cost  functions  using  the 
IDVD  (3rd)  method  are  displayed  in  Tables  23-25,  with  minimum  time  plotted  in  Figures 
123-124. 
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Figure  123.  Translational  discrepancies  for  propagated  IDVD  (3rd)  minimum  time 

solution. 
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Rotational  discrepancies  for  propagated  IDVD  (3rd)  minimum  time  solution. 
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Figure  124. 
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Table  23.  Value  of  terminal  conditions  at  the  final  time  for  IDVD  (3rd)  minimum  time  solution. 
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Table  24.  Value  of  terminal  conditions  at  the  final  time  for  IDVD  (3rd)  minimum  energy 

solution. 
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Table  25.  Value  of  tenninal  conditions  at  the  final  tune  for  IDVD  (3rd)  minimum  fuel  solution. 
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8.  Summary  of  Close  Approach  Results  for  Baseline  Maneuver 

Results  for  different  variations  of  the  IDVD  theme  pertaining  to  the  baseline 
maneuver  introduced  in  Chapter  III  and  discussed  throughout  this  dissertation  are 
summarized  in  Table  26. 


GPOPS 

IDVD  (Euler  Angles) 

(vary  endpoint  accelerations) 

Min 

cost 

Fuel 

Final  Time 

3.5086 

8.8943 

8.1 182 

3.7775 

10.0000 

9.0792 

Energy 

10.4471 

0.2445 

1.1587 

5.0869 

0.3495 

0.5738 

Fuel 

20.9419 

3.6941 

2.4863 

13.3627 

4.8462 

3.9329 

CPU  Time 

33,370.1 

86,040.7 

84,261.6 

10.6 

6.8 

11.0 

IDVD  (quaternions) 

(vary  endpoint  accelerations) 

IDVD  (quaternions) 

(vary  endpoint  jerk/accelerations) 

Min 
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Min 
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Fuel 
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Fuel 

Final  Time 

8.8868 

7.5181 

3.6582 

8.9137 

8.1442 

Energy 

0.2461 

0.3875 

6.2309 

0.2458 

0.4131 

Fuel 

13.4202 

3.7177 

2.8312 

15.0982 

3.7044 

2.7552 

CPU  Time  (sec) 

11.5 

16.1 

32.1 

86.6 

623.0 

405.5 

Table  26.  Summary  of  perfonnance  indices  and  computational  time  for  a  variety  of  trajectory 

generation  methods. 


For  the  cases  of  GPOPS  and  IDVD  employing  quaternions,  Figure  125  shows  the 
resulting  minimum  energy  control  histories  represented  in  a  single  plot.  Both  the  control 
histories  calculated  using  IDVD  match  well  with  the  solution  provided  by  GPOPS.  The 
similarities  of  the  solutions  are  further  reinforced  by  the  small  deviation  in  PI  with  the 
3rd  derivative  IDVD  solution,  IDVD  (3rd),  having  a  PI  within  0.65%  of  the  GPOPS 
solution  and  the  IDVD  (3rd-4th)  solution,  IDVD  (3rd-4th),  falling  within  0.54%. 


152 


Figure  125.  Control  history  of  different  methods  for  the  minimum  quadratic-control 

problem  formulation. 


The  overlaid  control  histories  for  the  minimum  time  case  are  shown  in  Figure  126. 


Figure  126.  Control  history  of  different  methods  for  the  minimum  time  problem 

formulation. 


153 


For  the  minimum  time  case,  the  reduced  performance  of  the  IDVD  solutions  is 
mainly  due  to  the  inability  to  implement  a  bang-bang  control  strategy.  This  is  because 
the  IDVD  results  in  a  smooth,  continuous  function  by  construction.  Although  the  IDVD 
(3rd-4th)  strategy  is  more  flexible,  having  the  added  capability  to  modify  the  endpoint 
jerk  values,  it  still  cannot  implement  the  discontinuous  switch  of  control  values  inherent 
to  a  minimum  time  solution  for  this  type  of  problem  and  shows  only  minor  improvement 
over  the  IDVD  (3rd)  strategy.  Although,  previous  work  calculating  bang-bang  solutions 
needed  to  generate  an  extra  control  of  layer  in  order  to  have  the  required  smooth  control 
history  for  interpolation  (Hurni  2009).  The  PI  of  the  IDVD  (3rd)  solution  was  within 
6.5%  of  the  GPOPS  solution  and  the  IDVD  (3rd-4th)  solution  was  within  4.3%. 

The  minimum  fuel  solutions  are  shown  in  Figure  127.  Although  the  IDVD 
solutions  have  a  slightly  higher  PI,  they  avoid  the  rapid  switching  characteristics 
portrayed  by  the  GPOPS  solution  inherent  to  the  bang-off-bang  control  nature  associated 
with  the  cost  function  described  by  Equations  (59)  and  (62).  The  solution  provided  by 
IDVD  is  also  smooth  through  a  region  that  GPOPS  provides  a  highly  oscillatory  solution 
incapable  of  being  implemented.  The  trade  off,  of  course,  is  the  slightly  higher  PI 
associated  with  the  solution,  13.9%  increase  for  IDVD  (3rd)  and  10.8%  for  IDVD  (3rd- 
4th). 
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Figure  127.  Control  history  of  different  methods  for  the  minimum  fuel  problem 

formulation. 
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VI.  FURTHER  SIMULATIONS  AND  ANALYSIS 


This  chapter  analyzes  the  rendezvous  final  approach  problem  in  several  different 
ways.  First,  the  methods  developed  in  this  dissertation  are  applied  to  previously 
published  examples.  Concepts  for  closed-loop,  real  time  operations  are  derived  and 
demonstrated  in  simulation.  Finally,  the  Inverse  Dynamics  in  the  Virtual  Domain  method 
is  applied  to  a  variety  of  scenarios  with  different  initial  conditions. 

A.  COMPARISON  WITH  PREVIOUS  STUDIES 

1.  Problem  Formulation  in  Two  Dimensions 

Simulations  were  run  to  compare  results  with  those  published  in  previous 
literature.  Specifically,  the  planar  case  of  Ma  (Ma  et  al.  2007)  that  is  based  on  the 
Sakawa  Shindo  (SS)  algorithm  (Sakawa  and  Shindo  1980).  This  article  attempts  to  solve 
the  planar  motion  minimum  time  problem  of  a  spacecraft  matching  position  and 
orientation  with  a  docking  point  that  is  rotating.  The  state  vector  is  defined  as: 

V  =  U  yr  0r  vxr  vyr  0)rf  (147) 

The  subscript  r  indicates  the  relative  motion  of  the  chaser  with  respect  to  the 
RSO,  observed  from  the  targets  body  fixed  frame.  There  are  two  force  controls  u\  and  ui, 
respectively  fixed  in  the  x  and  y  chaser  spacecraft  body  frame  and  one  control  torque  uj, 
fixed  along  the  z  axis.  The  controls  are  bounded  such  that:  -1  <  u  <  1 .  The  initial  and 
final  conditions  are  given  as: 


xr0  =  (10  10  n  12  1  1  1)T 
x;/  =  (1  0  0  0  0  0)T 


(148) 


Since  the  algorithm  needs  a  final  time,  t f ,  as  in  input,  the  method  is  to  solve  a 

minimum  fuel  problem,  reducing  t  f  at  each  iteration.  This  iteration  of  iterations  method 

can  be  computationally  expensive,  and  the  inability  to  find  a  solution  using  the  SS 
algorithm  does  not  guarantee  that  one  does  not  exist.  Furthennore,  relative  motion  due  to 
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the  rotating  reference  frame  (not  inertially  fixed)  and  the  mean  motion  due  to  Hill’s 
equation  are  neglected.  For  more  information  on  the  SS  algorithm,  the  reader  is  directed 
to  (Sakawa  and  Shindo  1980;  Sakawa  1999).  Although  the  algorithm  claims  to  satisfy 
Pontryagin’s  Minimum  Principle  at  every  iteration,  no  infonnation  about  the  costates  or 
transversality  conditions  are  supplied.  The  final  minimum  time  solution  arrived  at  by  this 
method  is  8.14  seconds  using  200  nodes  for  numerical  integration  by  the  Heun  method 
(Chyba,  Leonard  and  Sontag  2000).  The  computation  time  was  claimed  to  take  several 
minutes,  but  not  stated.  Instead,  the  author  stated  that  the  computational  time  was  not  the 
focus  of  the  research,  only  the  computation  of  a  feasible  solution. 

2.  Direct  Method  (GPOPS)  Formulation  and  Results 

The  same  problem  is  solved  by  restricting  the  GPOPS  and  IDVD  formulations 
derived  in  previous  sections  to  two  dimensions.  The  resulting  solution  for  GPOPS 
provided  a  minimum  time  for  the  maneuver  to  be  7.6695  seconds.  The  resulting  plot  of 
chaser  x  and  y  coordinates  fixed  are  presented  in  the  RSO  body  frame,  as  in  Ma  (2007),  is 
shown  in  Figure  128.  The  control  histories  are  shown  in  Figure  129,  as  well  as  the 
switching  conditions  developed  in  Chapter  III.B.  The  controls  obey  the  bang-bang  nature 
stated  in  Ma  (2007)  and  Chapter  III  of  this  dissertation,  which  are  dictated  by  the 
switching  functions.  The  controls  presented  in  Ma  are  said  to  approach  the  bang-bang 
structure  but  appear  rather  smooth.  Figures  130  and  131  show  the  time  history  of  the 
endpoint  conditions,  the  difference  in  position  and  velocity  of  the  docking  points,  as  well 
as  angular  rate  and  orientation  of  the  vehicles,  converging  to  zero.  The  transversality 
conditions  and  Hamiltonian  shown  in  Figures  132  and  133  further  reinforce  the  optimal 
nature  of  the  control  for  the  minimum  time  cost.  Still,  the  major  drawback  of  this 
formulation  and  method  is  that  the  calculation  time  of  the  solution  was  1,416.2  seconds 
(23.6  minutes)  for  only  a  50  node  solution  (not  shown),  and  the  calculated  solution  for 
200  nodes,  shown  in  Figures  128-133,  required  40,659.0  seconds  of  computational  time. 
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Figure  128.  2D  Case  (GPOPS):  optimal  minimum  time  trajectory  of  the  chaser  in  the  RSO 

body  frame. 
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Figure  129.  2D  Case  (GPOPS):  optimal  minimum  time  control  history. 
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Figure  130.  2D  Case  (GPOPS):  history  of  the  translational  endpoint  conditions. 
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Figure  131.  2D  Case  (GPOPS):  history  of  the  attitude  endpoint  conditions. 
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Figure  132.  2D  Case  (GPOPS):  history  of  the  transvervsality  conditions. 
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Figure  133.  2D  Case  (GPOPS):  history  of  the  Ffamiltonian. 


For  comparison  purposes,  the  GPOPS  method  was  run  with  a  minimum  fuel  cost 
based  on  Equation  (93)  but  with  a  multiplier  of  0.005  (Ma  et  al.  2007)  and  a  maximum 
final  time  of  8.14  seconds.  The  solution  took  1,953.1  seconds  (32.6  minutes)  to  calculate 
a  50  node  solution  and  provided  a  PI  of  0.0574  (compared  to  the  0.0788  PI  result  from 
the  SS  based  method). 
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3. 


IDVD  Formulation  and  Result 


The  IDVD  method  was  also  used  to  calculate  the  trajectory  based,  shown  in 
Figure  134,  on  the  2D  problem  posed  by  Ma  (Ma  et  al.  2007).  The  associated  control 
history  is  shown  in  Figure  135.  The  endpoint  conditions  are  shown  in  Figures  136  and 
137.  Again,  the  same  formulation  and  methodology  was  used  as  in  Chapter  V.A,  varying 
acceleration  at  the  endpoints,  but  the  problem  was  constrained  to  planar  motion.  The 
resulting  costs  and  plots  were  generated  using  200  points  (the  resulting  IDVD  is  analytic, 
so  any  number  of  nodes  or  waypoints  can  be  used  to  present  the  final  solution).  This 
final  time  calculated,  8.7525  seconds,  is  higher  than  both  the  previous  methods,  but  the 
solution  only  took  10.3  seconds  to  compute.  When  the  IDVD  (3rd-4th)  method  is  used, 
the  final  time  to  complete  the  maneuver  is  decreased  to  8.2086  seconds,  but  the 
computational  time  to  arrive  at  the  solution  increase  to  214.7  seconds. 


12 


IDVD  3rd 
IDVD  3rd-4th 


10 


2 

LL 


"O 

O 

CD 

O 

CO 

cc 

■2 

o 


0 


-2 


0 


2 


4  6  8  10  12  14 

X  coord  RSO  Body  Frame 


Figure  134.  2D  Case  (IDVD):  minimum  time  trajectory  of  the  chaser  in  the  RSO  body 

frame. 


160 


1 


CO 

3 


* 

* 

t 

/ 

-IDVD  3rd 
-IDVD  3rd-4th 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

- 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

\ 

% 

\ 

\ 

- 

N 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

Time  (s) 


Figure  135.  2D  Case  (IDVD  3rd):  minimum  time  control  history. 
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Figure  136.  2D  Case  (IDVD  3rd):  history  of  the  translational  endpoint  conditions. 
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Figure  137.  2D  Case  (IDVD  3rd):  history  of  the  attitude  endpoint  conditions. 

The  increased  flexibility  of  the  IDVD  (3rd-4th)  basis  functions  permits  the 
solution  to  behave  more  like  the  bang-bang  structure  of  the  optimal  solution,  allowing  for 
a  reduced  final  time  for  the  maneuver  to  complete. 

4.  Comparisons  of  Trajectory  Generation  Methods  for  2D  Example 

Table  27  shows  a  comparison  between  methods  for  the  minimum  time  solution. 
200  nodes  were  chosen  for  the  calculation  and  a  factor  of  0.005  was  multiplied  by  the 
minimum  fuel  cost  to  be  consistent  with  previous  literature  (Ma  2007).  Although  IDVD 
(3rd)  has  a  greater  final  maneuver  time,  the  computational  time  is  drastically  reduced. 
The  computational  time,  simply  stated  as  several  minutes,  and  final  minimum  energy  cost 
of  the  SS  solution  presented  by  Ma  (2007),  were  not  reported.  The  solution  providing  the 
best  minimum  time  cost  was  the  bang-bang  solution  solved  for  by  GPOPS,  but  the 
solution  came  at  an  extensive  computational  price. 
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GPOPS 

IDVD 

(3rd) 

IDVD 
(3rd -4th) 

Ma  (SS) 

cost 

Min 

Time 

Final  Time 

7.6695 

8.7525 

8.2086 

8.1400 

Energy 

11.3266 

5.9307 

6.9566 

N/A 

Fuel* 

0.1135 

0.0779 

0.0824 

0.0788 

CPU  Time  (sec) 

40,659.0 

10.3 

214.7 

N/A 

Table  27.  Summary  of  performance  indices  and  computational  time  for  a  2D  scenario.  *A  factor 
of  0.005  was  multiplied  by  the  minimum  fuel  cost  to  be  consistent  with  the  fonnulation 

presented  by  Ma  (2007). 

B.  CYCLICAL  NATURE  OF  PROBLEM  SOLUTION  PERFORMANCE 
INDEX 

The  performance  index  of  rendezvous  of  a  spacecraft  with  a  tumbling  object 
exhibits  a  cyclical  nature  as  seen  by  Figures  138  and  139  for  two  scenarios.  Figure  138 
shows  the  cyclic  nature  of  the  PI  with  respect  to  a  time  value  associated  with  waiting  to 
start  the  maneuver.  If  the  initial  angular  velocity  of  the  RSO  occurs  around  a  principal 
moment  of  inertia,  then  the  motion  of  the  docking  point  is  periodic,  containing  only 
circular  motion  in  the  plane  perpendicular  to  the  angular  velocity  vector.  A  simple 
example  of  fixing  A  =1  for  the  entire  maneuver  and  then  varying  the  wait  time  before 
maneuver  start  (in  essence  is  the  same  as  varying  the  initial  orientation  of  the  RSO), 
shows  that  the  PI  associated  with  completing  the  maneuver  is  periodic  with  time.  Figure 
139  shows  the  cyclic  cost  nature  for  the  case  of  the  docking  point  motion,  not  being 
constrained  to  a  plane  perpendicular  to  the  initial  angular  velocity  vector  of  the  RSO. 
Such  is  the  case  of  the  RSO  having  a  nonidentity  inertia  matrix  and  initial  angular 
velocity  vector  that  is  not  coincident  with  a  principal  moment  of  inertia,  the  PI  is  cyclic, 
but  does  not  repeat.  For  this  case,  the  time  value  shown  is  the  total  time  of  the  maneuver. 
For  both  cases,  there  is  an  inherent  cyclical  nature  to  the  PI  pertaining  either  the  wait  time 
until  commencing  the  maneuver  or  the  total  time  of  the  maneuver.  This  added 
complexity  demonstrates  the  increased  influence  of  the  maximum  final  time  on  the 
solution.  For  any  given  circumstance,  allowing  more  time  to  perform  the  maneuver  may 
not  necessarily  decrease  the  PI  for  minimum  fuel  or  minimum  energy  maneuvers  as 
would  be  the  case  if  rendezvous  was  to  a  fixed  point  in  space.  For  this  reason,  initial 

163 


guesses  should  be  chosen  judiciously  and  when  possible,  as  stated  in  Chapter  V.B,  and 
should  provide  a  feasible  solution.  This  way,  since  the  kinematic  and  dynamic  equations 
are  always  satisfied  though  IDVD,  there  is  always  a  solution  to  implement  at  any  given 
time  should  the  user  want  to  terminate  the  optimization  routine. 


Figure  138.  Energy  and  fuel  costs  for  rendezvous  with  a  tumbling  RSO  with  symmetric 
inertia  matrix.  The  time  shown  is  a  wait  time  until  the  start  of  the  maneuver. 


Figure  139.  Energy  and  fuel  costs  for  rendezvous  with  a  tumbling  RSO  with  an  asymmetric 

inertia  matrix.  The  time  shown  is  the  total  maneuver  time. 
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C.  FLEXIBILITY  AND  ROBUSTNESS  OF  IDVD  TECHNIQUE 

The  purpose  of  this  section  is  not  to  generate  an  exhaustive  list  of  potential 
advantages  and  implementations  of  the  IDVD  methods  described  in  the  previous  section, 
but  to  further  educate  the  reader  on  concepts  in  order  to  facilitate  new  ideas.  The 
minimum  energy  example  was  selected  for  further  analysis  and  implementation  based  on 
its  nature  to  minimize  control  effort  over  the  entire  maneuver.  It  has  inherent  fuel¬ 
minimizing  qualities  as  well  as  other  attractive  attributes.  The  minimum  energy  cost 
based  maneuver  does  not  have  the  trait  of  commanding  maximum  (or  saturated)  controls 
during  the  tenninal  phase  of  docking.  The  extra  control  margin  makes  it  the  most 
responsive  and  safe  if  sudden  changes  in  trajectory  are  needed,  a  potential  necessity  if 
closed  loop  implementation  is  to  be  considered.  Having  low  control  action  also  lowers 
potential  effects  of  plume  impingement  as  opposed  to  maximum  thrust.  These  combined 
attributes  make  it  the  preferred  choice  for  implementation. 

1.  Real  Time  Trajectory  Reshaping 

Because  the  IDVD  solution  is  analytic  in  nature  and  incorporates  prescription  of 
the  endpoint  values  as  well  as  their  derivatives,  it  allows  for  the  possibility  of 
implementation  in  a  closed-loop  fashion.  If,  for  example,  the  minimum  energy  trajectory 
(which  does  not  saturate  the  actuators  at  any  point)  solution  was  calculated  using  the 
IDVD  method,  the  resulting  varied  parameters  would  vary  the  higher  level  derivatives  at 
the  endpoints,  since  the  state  values  and  their  respective  first  derivatives  were  set  based 
on  the  projected  state  of  the  RSO  at  the  final  time  as  well  as  the  initial  conditions  of  the 
chaser  craft.  Even  after  a  solution  is  calculated,  if  updated  information  about  the  current 
states  of  the  RSO  were  obtained  and  the  final  state  values  of  the  RSO  were  to  change 
(due  to  unmodeled  dynamics  or  disturbances),  this  infonnation  can  be  reintroduced  into 
the  trajectory  generation  method  to  ensure  that  the  endpoints  of  the  trajectory  met  the 
desired  conditions.  Mathematically  speaking,  the  varied  parameters,  already  solved  for, 
would  not  change,  but,  keeping  the  conditions  at  the  beginning  of  the  trajectory  fixed,  the 
specified  conditions  at  the  ending  point,  tf,  would  be  tweaked  in  order  to  match  the  best 
projected  conditions  of  the  RSO  docking  point,  resulting  in  new  coefficients  for  the 
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polynomials  and  therefore  instantly  reshaping  the  trajectory.  The  new  trajectory  may  not 
be  optimized  with  respect  to  the  slight  changes  in  state  variables,  but  it  would  be  the  best 
solution  based  on  the  given  constraints  (having  a  feasible  trajectory  that  ends  with  the 
desired  conditions  on  the  chaser  and  its  docking  point)  without  going  through  another 
iteration  process.  This  is  in  contrast  to  current  iterative  real  time  closed  loop  optimal- 
control  studies  that  are  performed  in  simulation,  (McFarland  et  al.  2009),  and  assume 
solutions  can  be  obtained  in  whatever  specified  update  rate  is  needed,  stopping  the 
simulation  to  compute  the  optimal  control,  which  can  take  several  minutes  to  hours 
depending  on  computer  perfonnance  and  desired  solution  resolution.  The  following  two 
examples  are  used  to  highlight  the  noniterative  reshaping  trait  of  the  proposed  IDVD 
controller.  In  the  first  scenario,  the  inertia  matrix  of  the  RSO  is  thought  to  be  known  as  I 
=  diag([3, 1 ,2]);  therefore,  it  is  used  to  project  the  states  of  the  RSO  at  tf  based  on  the 
initial  conditions  of  the  RSO  and  Euler’s  equations  for  rotational  motion.  For  this 
exercise,  we  assume  the  solution  of  the  varied  parameters  that  optimize  the  rendezvous 
maneuver  is  already  known  and  the  conditions  are  such  that  the  problem  is  the  same  as 
discussed  in  previous  sections  based  on  Table  1.  The  only  difference  is  that  now  the 
actual  inertia  of  the  RSO  is  I  =  diag  ([1,2,3]).  The  second  scenario  involves  correct 
knowledge  of  the  RSO  inertia,  but  the  presence  of  an  unknown  constant  thrust  in  the 
body  frame  (such  as  a  thruster  that  is  failed  on).  New  (ideal)  RSO  state  infonnation  is 
available  at  a  rate  of  10  Hz,  allowing  the  trajectory  generator  to  reshape  in  real  time 
based  on  the  new  projected  endpoint  values.  A  block  diagram  of  the  process  is  shown  in 
Figure  140. 

In  the  previous  sections,  the  IDVD  Solver  coupled  with  the  Trajectory  Generator 
was  treated  as  the  same  system.  It  should  be  noted  that,  in  fact,  they  can  be  separated  in 
order  to  exploit  specific  advantages  to  the  IDVD  method.  The  key  difference  between 
the  IDVD  Solver  and  the  Trajectory  Generator  is  the  IDVD  Solver  iterates  on  the  varied 
parameters  (the  higher  level  derivatives  of  a  given  trajectory  at  the  endpoints)  and  uses 
the  Trajectory  Generator  to  solve  for  the  polynomial  coefficients  and  generate  the  spatial 
and  time  trajectories  of  the  system  states  (as  well  as  the  controls).  From  there,  a  PI  can 
be  associated  with  the  recently  calculated  trajectory  and  the  IDVD  Solver  can  iterate  on 


166 


those  varied  parameters  to  reduce  that  PI.  In  summary,  the  entire  IDVD  makes  calls  to 
the  Trajectory  Generator.  The  Trajectory  Generator,  when  examined  alone,  takes  the 
endpoint  conditions  (position  and  velocity)  and  higher  order  derivatives  (eg.  acceleration) 
for  attitude  and  translation  and  solves  for  the  polynomial  coefficients  using  Equations 
(12 1)— (128)  for  attitude  and  Equations  (1 36)— ( 139)  for  translation.  The  idea  being  that 
once  the  IDVD  Solver  has  values  for  the  varied  parameters,  the  most  current  information 
on  the  RSO  docking  point  can  be  used  to  reshape  the  trajectory,  ensuring  endpoint 
conditions  are  met.  Furthermore,  it  can  be  implemented  in  such  a  manner  that  also 
provides  the  controls  necessary  to  track  this  trajectory,  using  the  inner-loop  controller  to 
take  care  of  small  errors  (Yakimenko  et  al.  2008). 


Trajectory 
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Controller 


Initial 
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Updated  States 


Figure  140.  Block  diagram  of  real  time  trajectory  reshaping  concepts  using  IDVD. 

a.  Solution  Results  with  Inertia  Uncertainty 

The  trajectory  resulting  from  the  closed-loop  implementation  is  shown  in 
Figure  142.  For  this  case,  the  trajectory  generator  has  incorrect  information  about  the 
inertia  of  the  RSO,  specifically  it  is  using  a  diag  ([3,1,2])  inertia  matrix  (the  same  as 
previous  examples)  when  the  actual  inertia  matrix  is  diag  ([1,2,3]).  This  is  a  more 
extreme  case  of  misidentification  and  was  chosen  to  best  illustrate  the  concept.  In 
actuality,  the  inertia  matrix  information  of  the  RSO  employed  by  the  trajectory  generator 
would  most  likely  be  more  accurate,  or  simply  use  an  identity  inertia  matrix. 
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In  Figure  141,  the  actual  trajectory  of  the  RSO  docking  point  is  illustrated 
by  a  dotted  red  line  with  a  circle  marker  highlighting  a  subset  of  points  at  which  updates 
are  provided  to  the  RSO.  The  dash/dot  lines  illustrate  some  of  the  projected  trajectories 
of  the  RSO  docking  point  based  on  the  most  current  angular  position  and  velocity 
information,  propagated  with  the  (incorrect)  RSO  inertia  matrix  currently  used  by  the 
trajectory  generator.  The  dotted  line  shows  the  projected  endpoint  of  the  RSO  docking 
point  as  it  progresses  though  each  iteration.  The  projected  RSO  docking  point  and  the 
actual  RSO  docking  point  converge  to  the  same  value  as  at  tf,  since  the  amount  of  time 
for  the  trajectory  generator  to  project  into  the  future,  based  on  incorrect  inertia 
information,  where  the  RSO  docking  point  will  be  decreases.  The  thin  colored  lines 
show  a  subset  of  the  reshaped  trajectories  at  each  timestep.  The  blue  dots  show  the 
resulting  overall  trajectory  for  the  chaser  vehicle,  overlaid  with  intermittent  models 
showing  the  resulting  attitude  of  the  vehicle. 
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x  (radial)  (m) 


Figure  141 .  Results  of  the  final  trajectory  for  the  close  approach  example  having  a  misidentified  inertia  matrix  IDVD  real  time  reshaping. 
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Figure  142  shows  a  close-up,  expanded  view  of  the  previously  described  trajectories. 
Information  shown  is  the  actual  trajectory  (red,  circle),  the  projected  endpoint  (black, 
dash)  and  a  sample  of  the  overall  projected  trajectories  based  on  state  updates  (green, 
dash/dot).  Figure  143  illustrates  the  evolution  of  the  current  trajectory  over  time  and  the 
position  of  the  chaser  CM. 
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Figure  142.  Exploded  view  of  the  RSO  docking  point  information  for  a  misidentillcd 

inertia  matrix  example. 
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Figure  143.  Spatial  views  of  the  resulting  trajectories  for  the  close  approach  example 
having  a  misidentified  inertia  matrix  IDVD  real  time  reshaping. 

b.  Solution  Results  with  Unknown  Constant  Torque 

For  the  next  example,  it  is  assumed  the  inertia  matrix  of  the  RSO  is  known 
by  the  trajectory  generator,  but  an  unknown  torque  exists,  acting  on  the  vehicle  in  the  x- 
body  frame.  This  can  be  thought  of  as  a  thruster  that  is  stuck  in  the  On  mode,  thus 
generating  the  disturbance.  The  resulting  plots,  similar  to  those  in  the  previous  section 
are  shown  in  Figures  144  and  145.  Figure  144  shows  several  of  the  reshaped  trajectories, 
along  with  the  overall  trajectory  for  the  chaser  vehicle,  while  Figure  145  shows  the  time 
history  for  each  translational  component  of  the  Chaser  vehicle,  as  well  as  the  evolution  of 
the  reshaped  trajectories.  As  in  the  previous  example,  the  projected  and  actual  endpoint 
of  the  RSO  states  converge  to  zero  as  the  chaser  vehicles  converges  to  the  final  state. 
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Figure  144.  View  of  the  close  approach  trajectory  example  with  an  unidentified,  constant  torque  using  IDVD  real  time  reshaping. 
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Figure  145.  Spatial  views  of  the  resulting  trajectories  for  the  close  approach  example 
having  an  unidentified  constant  torque  using  IDVD  real  time  reshaping. 

2.  Resolving  and  Reoptimizing  the  Problem  in  Real  Time 

If  the  solution  of  the  IDVD  could  be  computed  in  real  time,  with  the  definition  of 
real  time  based  on  the  requirements  of  the  user/customer  and  the  estimation  attributes  of 
the  sensor  system,  the  problem  can  be  recomputed  as  the  trajectory  is  traversed.  An 
illustrative  example  of  this  concept  is  given  in  Figure  146.  In  this  case,  a  switch  is 
installed  that  allows  the  initial  conditions  to  be  updated  based  on  the  most  current  system 
states.  The  problem  can  then  be  resolved  using  IDVD  and  the  new  resulting  trajectory 
immediately  updated.  The  key  difference  between  this  implementation  and  the  reshaping 
approach  is  that  the  trajectory  is  resolved  and  reoptimized  based  on  the  most  current 
states  of  the  RSO  and  the  chaser  vehicle.  Even  though  for  this  case  the  trajectory  is 
resolved,  one  of  the  greatest  benefits  between  the  IDVD  method  is  still  that  if  the 
endpoint  conditions  change,  the  IDVD  can  use  this  information  to  tweak  the  trajectory. 
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Other  methods,  specifically  pseudospectral,  must  resolve  the  entire  problem  in  order  to 
obtain  a  trajectory  that,  when  the  controls  are  propagated,  would  finish  the  maneuver  in 
the  correct  position  with  the  desired  conditions. 


Figure  146.  Block  diagram  of  real  time  trajectory  recalculation  concepts  using  IDVD. 

This  concept  is  demonstrated  on  the  same  example  stated  in  Chapter  VII.D.l.b, 
employing  the  inertia  uncertainty.  The  trajectory  generator  is  hampered  by  not  knowing 
the  true  inertia  of  the  RSO,  but  is  receiving  state  updates  and  recalculating  the  trajectory 
at  a  rate  of  1Hz.  This  allows  the  current  trajectory  to  be  based  on  the  most  current  known 
states  of  the  chaser  and  the  RSO.  The  switch  allows  for  the  trajectory  generation  method 
to  use  the  rapid-reshaping  concept  described  in  Chapter  VII. D.l  with  a  10  Hz  update  rate 
for  the  final  1  second,  or  endgame  phase,  of  the  trajectory  since  uncertainties  in  guidance 
may  need  to  be  corrected  much  more  rapidly  (Vaidyanathan  et  al.  2001).  This  also 
demonstrates  the  IDVD’s  ability  to  be  integrated  with  other  guidance  systems  to 
construct  the  best  overall  system. 

Figure  147  shows  the  plots  of  each  recalculation  (top),  as  well  as  the  overall 
resulting  trajectory  (Bottom).  The  multicolored  segmented  lines  in  the  top  figure 
illustrate  each  trajectory  recalculation,  occurring  at  1 -second  intervals  shown  by  the  blue 
dots.  The  segmented  line  on  the  lm  sphere  shows  black  “x”  marks  representing  the 
current  projected  location  of  the  docking  point  at  the  calculated  tf  for  each  optimized 
trajectory. 
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Figure  147.  (Top)  Overlaid  recomputed  trajectories.  (Bottom)  Resulting  overall  trajectory 

based  on  the  most  current  recalculation. 


3.  Specific  Waypoint  Spacing  in  the  Time  Domain 

Another  advantage  to  the  specific  IDVD  formulation,  and  an  enabling  attribute 
that  allowed  the  implementation  of  the  closed-loop  implementation  concepts  from  the 
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previous  sections,  is  the  freedom  to  choose  node  or  waypoints.  In  previous  work 
(Bevilacqua,  Romano  and  Yakimenko  2009),  integrating  the  inverse  of  speed  factor  that 
was  based  on  switching  times  led  to  the  user  to  be  at  the  mercy  of  node  spacing  that  is 
dependent  upon  the  history  of  the  speed  factor  itself.  This  could  lead  to  gaps  in  control 
and  trajectory  information  that  make  the  solution  incapable  of  implementation.  This  is 
also  seen  in  application  of  pseudospectral  direct  methods  (Humi  2009)  while  using  a 
commercial  program  for  direct  optimization  of  trajectories  and  having  to  add  layer  of 
control  (and  complexity)  in  order  to  correctly  interpolate  control  information  in  between 
the  node  points  of  the  solution.  Because  of  this,  for  the  problem  studied  concerning  a 
wheeled  rover-like  vehicle,  the  user  is  forced  to  use  velocity  and  heading  as  a  control 
vector,  a  layer  above  the  intuitive  controls  of  acceleration  and  angular  rate.  An  example 
where  the  constraints  on  the  previous  problem  are  relaxed  to  illustrate  how  the  nodes  of 
the  solution  get  bunched  together  in  regions  where  the  speed  factor  is  the  greatest  is 
considered.  This  is  because  the  node  spacing  is  uniform  with  respect  to  the  virtual 
argument  t,  as  shown  in  the  bottom  plot  of  Figure  148. 
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Figure  148.  Control  history  of  a  representative  trajectory  with  uniform  node  spacing  in  the 

virtual  domain. 
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The  question  posed  is  what  dispersion  of  node  points  in  the  virtual  domain  will 
guarantee  uniform  node  spacing  in  the  time  domain.  For  the  methods  described  in 
Chapters  III  and  IV,  the  speed  factor  was  chosen  such  it  has  an  analytic  expression  for  the 
integral  of  its  inverse,  which  maps  a  specific  value  of  t  to  the  time  domain  as  shown 
below: 


t  = 


(149) 


This  results  in  the  ability  to  express  z  as  a  function  of  t  and  vice  versa. 

Once  the  final  time  is  obtained  from  the  initial  solution,  a  set  of  evenly  space 
nodes  (or  arbitrarily  spaced  nodes  based  on  any  inputs  from  the  user)  set  in  the  time 
domain  can  be  instantaneously  converted  to  the  virtual  domain.  This  new  set  of  nodes  in 
the  virtual  domain,  when  implemented,  leads  to  the  evenly  spaced  nodes  in  the  time 
domain  as  shown  in  Figure  149. 
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Figure  149.  Control  history  of  a  representative  trajectory  with  uniform  node  spacing  in  the 

time  domain. 
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D.  EXTENDED  EXAMPLES  USING  INVERSE  DYNAMICS  IN  THE 
VIRTUAL  DOMAIN 

Up  to  this  point,  the  same  initial  conditions  have  been  used  to  analyze  the  optimal 
solution  to  a  close  approach  problem,  compare  with  a  suboptimal  rapid-trajectory 
generation  technique,  and  demonstrate  potential  real  time  close-loop  applications.  The 
IDVD  (3rd)  and  IDVD  (3rd-4th)  method  generated  in  the  previous  chapter  is  applied  to 
several  pseudorandom  initial  chaser  conditions  and  scenarios.  First,  10  samples  were 
developed  to  apply  the  IDVD  technique.  From  these  10  samples,  the  initial  velocities  of 
the  spacecraft  were  assumed  to  be  zero  and  the  initial  orientation  of  the  chaser  is  assumed 
to  be  coincident  with  the  orbit  frame;  both  assumptions  were  taken  from  initial  conditions 
used  by  McCamish  (2007).  The  pseudorandom  initial  positions  of  the  chaser  are  chosen 
by  first  detennining  a  vector  them  multiplying  it  by  a  uniformly  distributed  random 
number  between  5  m  and  10  m  as  shown  in  Equation  (150). 

x  =  -  randx2  )cos  (rand  e)r and  r 

v  =  ^j(  1  -  rand^  j  sin  (randd  )randr  (150) 

z  =  rand /and, r 

with  rand  x  e  [0  1]  ,rand0e[ 0  In'),  randr  e  [5  10]  having  unifonn  random 

distributions.  The  desire  was  to  have  an  equally  distributed  set  of  points  in  range  (from 
5-10  m)  and  all  directions,  understanding  that  based  on  the  given  formula,  the  points 
would  not  necessarily  be  uniformly  distributed  about  the  volume  housing  those  points. 
The  solutions  to  the  first  ten  generated  initial  conditions  using  IDVD  (3rd)  are  illustrated 
in  Figure  150.  Again,  the  initial  conditions  of  the  RSO  were  kept  consistent  with  the 
previous  examples  from  this  manuscript,  having  a  I  =  diag[3  1  2]  and  initial  angular 
velocity  of  0.25  radians/second  in  the  y  and  z  body  frame  and  initial  q  =  [0  0  0  1].  For 
the  first  set  of  ten  shown  in  Figure  150,  the  average  computational  time  was  21.0 
seconds. 
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Figure  150.  Trajectories  for  sample  rendezvous  problem  for  the  first  10  random  starting 

chaser  positions  with  zero  initial  relative  velocity. 

Next,  over  1000  pseudorandom  initial  conditions  were  generated  with  varying 
initial  relative  velocity  as  well  as  position.  The  maximum  magnitude  of  the  initial 
velocity  was  limited  to  1  m/s  to  be  consistent  with  “relative  speed  limit”  constraints 
placed  during  close  proximity  operation  again  in  the  complementary  work  of  McCamish 
(2007)  as  well  as  having  the  body  frame  aligned  with  the  orbital  frame.  This  would  be 
consistent  with  a  GNC  concept  utilizing  McCamish’ s  Artificial  Potential  Function 
modified  algorithm  up  until  it  is  time  to  perform  the  close  approach  for  docking 
maneuver  when  the  RSO  has  nonzero  angular  rates.  The  solutions  to  the  first  10  samples 
are  shown  in  Figure  151.  Notice  the  presence  of  initial  velocity  for  the  chaser  vehicles 
results  in  a  modified  path  to  the  final  conditions.  This  is  due  to  the  fact  the  resulting 
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polynomial  trajectory  must  shape  itself  to  have  a  directional  component  at  the  beginning 
of  the  maneuver  that  will  accommodate  the  different  initial  conditions. 


S'  ? 


Figure  151.  Trajectories  for  sample  rendezvous  problem  for  the  first  10  random  starting 

chaser  positions  with  nonzero  initial  relative  velocity. 

The  complete  distribution  of  initial  conditions  is  shown  in  Figures  152-155.  The 
maximum  allowable  final  time  was  also  increased  to  15  seconds  for  the  ensuing 
calculations. 
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Figure  152. 


Distribution  of  the  1000  pseudorandom  initial  conditions  tested. 


in-track  (m)  in-track  (m)  cross-track  (m) 

Figure  153.  Planar  view  of  the  1 000  pseudorandom  initial  conditions  tested. 
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Magnitude  of  Initial  Range  from  RSO  (m) 


Figure  154.  Distribution  of  the  initial  chaser  range  for  the  1000  pseudorandom  samples 

generated. 
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Figure  155.  Distribution  of  the  initial  chaser  velocity  for  the  1000  pseudorandom  samples 

generated. 
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The  IDVD  (3rd)  method  was  run  using  the  1000  sample  and  the  mean  time  to  calculate 
the  solution  was  34.3  seconds.  The  mean  time  to  complete  the  maneuver  was  14.3081 
seconds  and  the  mean  energy  PI  was  0.2289.  The  same  1000  initial  condition  set  was 
also  run  using  the  IDVD  (3rd-4th)  method.  The  mean  of  the  energy  PI  was  0.21 12,  but 
the  average  run  time  was  228.9097  seconds.  A  summary  of  the  results  is  shown  in  Table 
28  for  initial  conditions  on  relative  position  and  velocity  pertaining  to  Figures  154  and 
155  and  Table  29  pertaining  to  initial  relative  positions  based  on  Figure  154  but  starting 
with  zero  initial  relative  velocity.  Table  30  shows  the  results  for  varying  the  initial 
conditions  on  position  and  velocity  of  the  chaser  vehicle  according  to  Figures  154  and 
155  but  also  varying  the  initial  angular  velocity  of  the  RSO  based  on  Figure  156. 
Although  the  IDVD  (3rd-4th)  method  results  in  an  overall  average  cost  over  the  sample 
of  initial  conditions  tested,  the  computational  time  was  substantially  more. 


IDVD 

(3rd) 

IDVD 

(3rd-4th) 

mean  values 

Min 

Energy 

Min 

Energy 

Final  Time 

14.3081 

14.5191 

Energy 

0.2289 

0.2112 

Fuel 

4.3124 

4.2516 

CPU  Time  (sec) 

34.3 

228.9 

Table  28.  Summary  of  perfonnance  indices  and  computational  time  for  IDVD  (3rd)  and  IDVD 
(3rd-4th)  method  using  1000  pseudorandom  samples  for  varying  the  initial  conditions 

of  the  chaser  position  and  velocity. 


IDVD 

(3rd) 

IDVD 

(3rd-4th) 

mean  values 

Min 

Energy 

Min 

Energy 

Final  Time 

14.3094 

14.6107 

Energy 

0.1631 

Fuel 

3.8544 

CPU  Time  (sec) 

32.9 

191.0 

Table  29.  Summary  of  perfonnance  indices  and  computational  time  for  IDVD  (3rd)  and  IDVD 
(3rd-4th)  method  using  1000  pseudorandom  samples  for  varying  the  initial  conditions 
of  the  chaser  position  but  having  the  initial  relative  velocity  be  equal  to  zero. 


183 


Magnitude  of  Initial  RSO  Angular  Velocity  (rad/s) 


Figure  156.  Distribution  of  the  initial  RSO  angular  velocity  for  the  1000  pseudorandom 

samples  generated. 


IDVD 

(3rd) 

IDVD 

(3rd-4th) 

mean  values 

Min 

Energy 

Min 

Energy 

Final  Time 

14.9716 

14.9799 

Energy 

0.2138 

0.1689854 

Fuel 

3.3867 

3.2472 

CPU  Time  (sec) 

20.0 

228.5 

Table  30.  Summary  of  performance  indices  and  computational  time  for  IDVD  (3rd)  and  IDVD 

(3rd-4th)  method  using  1000  pseudorandom  samples  for  varying  the  initial  conditions 
for  position  and  velocity  of  the  chaser  and  initial  angular  velocity  of  the  RSO. 
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VII.  CONCLUSION 


This  section  summarizes  the  contributions  and  findings  of  this  dissertation  and  is 
organized  as  follows:  First,  a  summary  of  the  research  is  presented.  This  is  followed  by  a 
summary  of  the  contributions  of  this  research  followed  by  a  discussion  of  the  advantages 
and  disadvantages  pertaining  to  the  different  methods.  Finally,  a  recommendation  for 
future  study  is  presented. 

A.  SUMMARY  OF  RESEARCH 

A  six  DoF,  20-state  model  of  two  spacecraft  rendezvous  is  developed,  one  of 
which  was  controlled,  the  other  considered  to  be  passively  tumbling.  A  solution  is 
obtained  for  the  problem  of  close  approach,  up  to  the  point  of  contact,  using  a  direct 
optimal  control  method.  The  solution  is  then  verified  as  optimal  by  way  of  an  indirect 
method  based  on  the  MP.  Next,  a  trajectory  generation  method  for  spacecraft 
reorientation  is  developed,  based  on  a  quaternion  construction  of  IDVD.  This  new 
construction  enables  implementation  of  an  IDVD  trajectory  generation  method  for  the 
problem  of  a  spacecraft  performing  a  close  approach  maneuver  to  a  tumbling  object. 
Finally,  the  advantages  of  the  IDVD  method  were  exploited  and  demonstrated  in 
simulated  scenarios  that  employ  closed-loop  feedback. 

B.  SUMMARY  OF  CONTRIBUTIONS 

Contributions  of  this  dissertation  are  considered  to  have  two  distinct  categories  of 
impact,  theoretical  contributions,  enhancing  the  overall  understanding  contribute  to  the 
overall  body  of  knowledge,  and  practical  contributions,  advancing  the  current  state  of  the 
art  and  moving  the  technology  closer  to  implementation.  They  are  both  summarized  in 
the  next  several  paragraphs. 

Summarizing,  the  minimum  quadratic-control,  minimum  fuel  and  minimum  time 
continuous  optimal  control  3D  spacecraft  rendezvous-docking  problems  are  fonnulated 
for  the  first  time  in  literature,  and  addressed  using  a  direct  collocation  method,  GPOPS. 
For  both  problems,  the  desired  optimal  trajectory  of  chaser  spacecraft  with  respect  to  a 
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tumbling  RSO  is  sought  such  that  the  desired  docking  points  match  in  position  and 
velocity.  Moreover,  the  solutions  obtained  were  verified  in  the  numerical  simulations 
based  on  the  Minimum  Principle.  Those  solutions  included  derivation  of  the  adjoint 
equations,  formulation  of  the  necessary  conditions  for  the  optimal  solution,  and  synthesis 
of  the  optimal  controls.  The  results  obtained  from  using  the  direct  collocation  method — 
in  this  case  the  open  source  software  of  GPOPS — are  very  close  to  those  obtained  using 
with  the  Minimum  Principle.  It  was  also  found  that  path  constraints  are  necessary  when 
solving  for  the  optimal  trajectory  in  order  to  prevent  undesired  collision  of  spacecraft. 
This  was  shown  by  the  active  path  constraint  that  results  in  discontinuous  costates  upon 
contact  with  the  constraint  boundary. 

As  expected,  the  GPOPS  direct  collocation  (pseudospectral)  method  was  found  to 
be  quite  reliable  and  yielded  the  results  relatively  fast.  Moreover,  these  results  were  also 
validated  via  propagation  of  the  states  based  on  the  system  dynamics,  and  optimal 
controls  were  found.  That  is  what  would  be  needed  for  possible  onboard  implementation 
of  the  developed  algorithms.  However,  the  GPOPS  direct  collocation  method  proved  to 
be  unable  to  produce  the  results  fast  enough  so  that  it  could  be  used  in  real  time.  Hence, 
in  the  nearest  future,  it  only  leaves  off-line  open-loop  option  to  be  possibly  used  on  the 
real  spacecraft.  But  even  then,  it  was  found  that  in  case  of  a  singular  control,  the  method 
returns  somewhat  infeasible  results  that  cannot  be  implemented.  The  main  simplifying 
hypothesis,  considered  in  the  numerical  simulation  presented  in  this  section  of  the 
research,  was  the  spherical  inertial  symmetry  of  both  chaser  and  target  spacecraft,  as  well 
as  thrust  control  limits  being  applied  in  the  orbital  frame  regardless  of  chaser  spacecraft 
orientation.  In  order  to  achieve  feasible  controls  in  real  time,  this  assumption  needs  to  be 
removed,  along  with  investigating  the  usage  of  different  combined  Pis  and  other  direct 
methods  based  on  inversing  of  the  dynamics  of  the  problem. 

This  research  then  continued  work  on  spacecraft  rendezvous  by  modifying  the 
previous  problem  formulation  to  include  body  mounted  translational  thrusters  on  the 
chaser  spacecraft,  as  well  as  extensions  on  the  final  conditions.  These  new  final 
conditions  involved  matching  the  velocity  and  position  of  the  spacecraft  docking  points, 
as  well  as  matching  angular  velocity  and  orientation.  Inertia  parameters  on  the  RSO  were 
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also  set  to  different  values,  removing  the  assumption  of  an  identity  inertia  matrix,  in  order 
to  stimulate  dynamic  angular  rates  in  all  three  axes.  The  problem  was  addressed  again 
using  the  most  plausible  pseudospectral  methods.  In  addition,  the  optimal  control 
structure  was  synthesized  and  analytical  expressions  for  the  transversality  conditions 
were  derived  (exploiting  the  Minimum  Principle)  to  compare  those  with  the  calculated 
solution.  In  this  section  of  research,  three  Pis  were  examined.  It  was  found  that  the 
minimum  time  solution  exhibits  the  expected  bang-bang  control  and  the  minimum  fuel 
solution  exhibits  bang-off-bang  characteristics.  Assuming  controls  in  the  body  frame  as 
opposed  to  the  orbital  frame  (as  in  the  previous  work  by  the  authors),  while  making  the 
problem  more  complex,  seems  to  alleviate  the  singular  control  that  appeared  in  the  z- 
translation,  and  therefore,  made  the  GPOPS  solution  feasible.  Another  finding  was  that 
instead  of  addressing  the  minimum  fuel  problem,  it  was  wiser  to  consider  a  minimum- 
energy  problem  instead.  The  reason  for  this  is  that  although  the  minimum  fuel  solution 
happened  to  be  about  30%  more  effective,  it  would  be  difficult  to  implement  the  control 
in  practice  because  of  the  bang-off-bang  control  structure.  The  minimum  energy  solution, 
which  lasts  longer  (but  follows  a  similar  path),  assures  smooth,  easy-to-follow  control 
time  histories.  Furthermore,  the  minimum  energy  solution  does  not  involve  having  the 
actuators  saturated  at  the  final  time  (a  characteristic  of  bang-off-bang  maneuvers  for 
minimum  fuel  or  bang-bang  for  minimum  time),  increasing  the  safety  of  the  maneuver  by 
having  control  margin  available  if  a  rapid  corrective  maneuver  is  needed  and  reducing  the 
likelihood  for  or  large  disturbances  due  to  plume  impingement.  The  minimum  time 
solution  (following  a  different  path)  does  not  contact  any  state  constraints  and  therefore 
has  continuous  costates,  while  the  minimum-control  and  minimum  energy  solutions 
contact  the  keep-out  zone  of  the  RSO,  while  approaching  for  docking.  Finally,  it  was 
shown  that  the  required  CPU  time  greatly  exceeds  that  of  the  maneuver  itself,  which  as  of 
now  excludes  pseudospectral  direct  collocation  methods  from  being  a  candidate  for 
onboard  application.  Instead,  exploitation  of  other  approaches,  willing  to  sacrifice  a 
fraction  of  the  PI  but  decrease  the  required  CPU  time  by  several  orders  of  magnitude, 
were  investigated,  as  summarized  below. 
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In  the  case  of  reorientation  maneuvers,  the  IDVD  method  with  the  novel 
quaternion  approximation  functions  proposed  in  this  dissertation  allows  computing 
feasible  solutions  fast  enough  to  be  used  onboard  satellites  for  on-line  computation  of 
slew  maneuvers.  Moreover,  because  of  the  smooth  controls'  histories  it  can  be 
implemented  in  the  control  schemes  involving  a  feed-forward  loop.  For  this  specific 
implementation  a  form  of  the  speed  factor  was  chosen  based  on  a  combination  of 
nonlinear  equations  as  opposed  to  integrating  dynamic  equations  based  on  having 
switching  times  of  the  controls  act  as  varied  parameters  (Yakimenko  2000).  Because  of 
this,  compared  to  the  true  time-optimal  solutions,  the  IDVD  trajectories  do  not  have  the 
capability  to  generate  a  bang-bang  control  solution,  which  results  in  a  slightly  worse  PI. 
However,  smooth  controls  benefit  other  mission  preferences  of  having  desired  rates  at  the 
endpoints  as  well  as  guaranteeing  the  endpoint  constraints  of  the  maneuver  are  always 
satisfied.  In  addition,  it  is  a  definite  advantage  in  rapidly  changing  acquisition  or  tracking 
scenarios  and  when  the  slewing  spacecraft  possesses  low  frequency  flexible  modes.  This 
formulation  is  later  applied  to  situations  where  the  attitude  is  coupled  with  other 
dynamics  such  as  translational  motion  in  rendezvous  and  docking  applications.  In  this 
case,  a  simple  eigenaxis  slew  would  not  meet  mission  criteria  such  as  matching  rotational 
motion. 

IDVD  is  further  applied  to  the  maneuver  that  has  both  translational  and  attitude 
motion,  such  as  the  case  with  the  focus  of  this  research  involving  close  approach  to  a 
tumbling  object.  The  IDVD  method  developed  is  highly  flexible,  allowing  for  matching 
any  mixed  set  of  endpoint  conditions,  generating  user  specified  waypoint  spacing,  and 
providing  a  mechanism  for  real  time  closed-loop  trajectory  reshaping.  Based  on  the 
trajectory  generation  results  for  the  benchmark  3D  scenario,  the  minimum-control 
solution  was  chosen  for  further  study  and  implementation.  This  was  due  to  its  inherent 
characteristic  to  conserve  fuel,  while  providing  specific  safety  qualities,  such  as  not 
commanding  actuator  saturation  (which  would  limit  its  ability  to  do  rapid  reactive 
maneuvers  in  case  of  closed  loop  implementation  of  higher-level  emergency  abort). 
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Within  the  IDVD  method,  there  are  still  variations  in  trajectory  parameterization 
techniques.  The  IDVD  3rd  method  employs  varying  the  3rd  derivative  at  the  endpoints 
of  the  trajectory  while  the  IDVD  3rd-4th  method  uses  both  the  3rd  and  4th  derivatives  as 
varied  parameters.  Because  of  this,  the  IDVD  3rd-4th  method  can  provide  a  final 
trajectory  that  is  lower  in  PI  but  at  the  cost  of  computational  time.  The  IDVD  3rd  method 
can  obtain  a  solution  that  is  within  10%  of  the  IDVD  3rd-4th  method  while  only  using 
15%  of  the  computational  time,  based  on  results  from  1000  samples  with  random  initial 
conditions. 

C.  CONSIDERATIONS  AND  RECOMMENDATIONS 

Trajectory  generation  by  direct  collocation  methods  have  the  benefit  of  providing 
a  mathematically  verifiable  optimal  solution.  The  major  drawback  is  the  excessive 
computational  time  needed  on  current  computing  systems.  The  resolution  of  the  solution, 
and  computational  time,  is  highly  dependent  on  the  amount  of  nodes  used  to  solve  the 
problem.  Furthermore,  the  user  cannot  specify  where  the  nodes  are  placed  and  this  leads 
to  complex  interpolation  schemes  or  adding  an  extra  layer  of  control,  increasing  the 
complexity  of  the  problem  formulation  in  order  to  have  sufficient  control  infonnation  to 
implement.  The  benefits  seen  by  IDVD  methods  are  the  rapid  computational  time  that 
allows  for  a  feasible  solution  to  be  generated,  potentially  onboard  a  spacecraft  and  in 
closed-loop.  Because  of  the  analytic  nature  of  its  solution,  high  node  trajectories  can  be 
generated  without  sacrificing  computational  time.  The  ability  to  reshape  itself  when 
provided  updated  information  on  the  RSO  states  and  the  docking  point  reinforce  the 
overall  robustness  of  the  IDVD  method  for  trajectory  planning.  By  increasing  the  order 
of  the  basis  polynomials,  thereby,  increasing  the  varied  parameters,  solutions  can  be 
obtained  that  better  minimize  the  PI,  but  at  the  cost  of  higher  computational  times. 
Furthermore,  the  node  spacing  of  the  control  solution  can  be  completely  specified  by  the 
user,  therefore,  eliminating  the  previously  discussed  error  incurred  by  interpolation  of 
potentially  complex  control  profiles  or  adding  extra  layers  of  control  because  of  the  lack 
of  flexibility  in  the  node  spacing  in  the  control  solution  set  by  the  direct  collocation 
methods,  which  are  not  guaranteed  to  be  continuous.  It  should  be  noted  that  using  the 

IDVD  formulation,  constraints  can  be  set  on  higher  derivatives  than  the  level  of  the 
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control  vector.  For  example,  constraints  on  the  jerk  profiles  can  be  enforced,  while  still 
using  the  acceleration  profile  as  the  control  vector.  This  attribute,  which  is  in  complete 
contrast  other  direct  methods,  is  based  on  the  analytical  formulation  of  the  trajectory.  For 
example,  to  set  constraints  on  jerk  using  direct  collocation  method  techniques,  the  control 
vector  would  have  to  at  least  be  based  on  jerk  (or  even  higher  derivative)  and  not 
acceleration,  again  adding  additional  levels  of  complexity. 

Drawbacks  to  the  IDVD  method  are  that,  while  they  are  optimized  to  give  the  best 
trajectory  based  on  the  constraints  of  the  basis  polynomials,  they  cannot  give  the  truly 
optimal  solution  described  by  the  infinite  dimensional  optimal  control  problem,  nor 
would  an  optimal  solution  be  verifiable  through  Pontryagin’s  MP  due  to  the  lack  of 
costate  information.  Also,  maneuver  duration  is  an  issue  that  currently  plagues  both  the 
direct  and  IDVD  methods.  For  IDVD  this  is  due  to  the  additional  coupling  of  the  attitude 
and  translational  motion  through  the  speed  factor.  If  the  speed  factor  is  made  small  to 
accommodate  constraints  on  translational  motion  and  there  is  a  nonzero  requirement  on  a 
higher  order  derivative  at  an  attitude  maneuver  endpoint  (note  this  issue  does  not  exist  for 
translational  maneuvers  or  rest-to-rest  maneuvers),  then  the  attitude  maneuver  may 
"wrap"  or  rotate  beyond  2n  in  order  to  perfonn  the  maneuver  to  meet  the  endpoint 
requirements.  This  is  because,  for  the  given  set  higher  order  derivative  in  the  time 
domain,  a  reduction  of  speed  factor  at  the  end  point  (which  may  be  needed  to 
accommodate  maneuvers  with  excessive  times)  would  lead  to  a  larger  value  for  the 
respective  higher  order  derivative  in  the  virtual  time  domain  since  X  appears  in  the 
denominator  of  the  conversion  multiplier  from  time  to  virtual  domains.  This  also 
translates  into  a  performing  a  larger  maneuver  in  quaternion  space  to  accommodate  the 
higher  derivative  value  in  the  virtual  domain.  While  the  maneuver  would  still  be 
completely  feasible,  the  extra  rotation  may  be  undesired  in  the  final  approach.  Even 
though  this  is  not  an  issue  for  the  implementation  discussed  in  this  dissertation,  it  would 
need  to  be  addressed  for  application  on  a  subset  of  missions  that  have  drastically  different 
timescales  associated  with  the  attitude  and  translational  maneuver.  Although  this  is 
easily  rectified  by  simply  setting  constraints  on  the  speed  factor  so  the  extra  rotation  does 
not  happen,  it  may  hinder  the  resulting  PI.  Another  approach  is  employing  a  trajectory 
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generation  scheme  for  maneuvers  starting  farther  away  (far  approach),  which  simply 
holds  attitude  fixed  in  the  orbital  frame  (McCamish  2007)  or  pointing  to  the  RSO,  then 
switching  to  the  final  approach  scheme  discussed  in  this  dissertation.  This  is  in  line  with 
operational  methodology  that  would  want  to  be  pointing  at  the  target  as  long  as  possible 
in  order  to  get  the  most  up-to-date  information  on  the  tumbling  RSO  before  performing 
the  rapid  final  approach  maneuver  at  the  last  possible  instant.  On  the  other  hand,  for  the 
direct  collocation  methods,  maneuvers  that  are  particularly  long  in  length  are  subject  to 
large  time  gaps  in  between  nodes,  where  no  state  or  control  infonnation  is  available. 
Extra  caution  needs  to  be  taken  when  employing  such  a  solution. 

Summarizing,  two  novel  guidance  strategies  for  final  approach  to  a  tumbling 
object  are  developed  in  this  dissertation.  One  is  based  on  Inverse  Dynamics  in  the 
Virtual  Domain,  employing  a  newly  fonnulated  quaternion  approach,  and  the  other  is  a 
novel  formulation  that  employs  direct  collocation  methods  for  optimal  trajectory 
generation.  While  the  direct  collocation  method  formulation  is  too  computationally 
expensive  to  be  employed,  it  provides  a  mathematically  verifiable  optimal  baseline  for 
the  maneuvers.  While  the  method  based  on  1DVD  cannot  exactly  match  the  optimal 
solution,  it  can  provide  a  feasible  that  is  optimized  stemming  from  set  of  basis  functions 
with  a  fraction  of  the  computational  time  of  the  direct  method.  This  and  other  attractive 
features  discussed,  allows  for  exploitation  of  several  real  time  implementation  concepts, 
ft  is  recommended  that  the  1DVD  method  be  employed  for  situations  where  the 
computational  cost  of  the  solution  outweighs  the  PI  associated  with  the  state  and  control 
histories.  The  direct  collocation,  or  pseudospectral,  method  should  be  used  in 
circumstances  where  computational  time  is  not  a  concern,  but  a  truly  optimal  solution 
based  on  the  state  and  control  variables  is  desired,  such  as  base  lining  scenarios  and 
finding  limits  of  specific  technologies. 
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D.  POSSIBLE  FUTURE  DEVELOPMENTS 

The  following  research  problems  remain  open  for  possible  follow-up  efforts: 

1.  Investigate  effects  of  different  sets  of  basis  function  for  trajectory 
generation  and  overall  effect  of  using  IDVD  methods  to  seed  initial 
guesses  for  pseudospectral  optimal  control  solvers. 

2.  Analyze  performance  of  trajectory  tracking  in  simulation. 

3.  Experimental  validation  on  real  autonomous  vehicles. 

4.  Integration  with  current  research  involving  navigation  and  target 
identification. 
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APPENDIX.  QUATERNION  PROPERTIES 


A  quaternion  is  a  four-dimensional  vector  that  is  used  to  express  orientation  of  a 
rigid  body.  It  is  composed  of  a  scalar,  5,  and  a  vector  portion,  v.  There  are  several 
representations  and  formulation  of  quaternions  and  their  associated  properties  in 
literature.  Throughout  this  dissertation,  the  properties  in  this  Appendix  will  be  employed 
(Titterton  and  Weston  1997;  Sarkka  2007). 


Quaternion  Expression: 


UJ 


(151) 


Quaternion  Length: 


(152) 


Quaternion  Logarithm: 


(153) 


Quaternion  Exponential: 


exp(q)  =  exp(s) 


(154) 
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Quaternion  Conjugate: 


/ 

-q2 

-q3 
V~q*J 


Quaternion  Inverse: 


Quaternion  Product: 


'  Pi' 
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Pi 
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q  *p  = 
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